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Abstract  of  Dissertation  Presented  to  the  Graduate  School 
of  the  University  of  Florida  in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of  Doctor  of  Philosophy 

FINITE  ELEMENT  ANALYSIS  OF  VISCOELASTICALLY 
DAMPED  COMPOSITE  STRUCTURES 

By 

Vidish  S.  Rao 
August  1991 

Chairman:  Dr.  Chang-Tsan  Sun 
Cochairman:  Dr.  Bhavani  V.  Sankar 

Major  Department:  Aerospace  Engineering,  Mechanics  and 

Engineering  Science 

The  main  objective  of  this  study  was  to  examine  the 
potential  for  using  passive  viscoelastic  damping  treatments  to 
improve  the  response  characteristics  of  composite  structures. 
An  efficient  finite  element  analysis  model  to  estimate  the 
response  characteristics  of  viscoelastically  damped, 
prestressed  laminated  composite  structures  is  presented. 
Specialized  structural  and  continuum  finite  elements,  ideally 
suited  to  analyze  such  structures,  were  formulated  and 
implemented  in  a stand-alone  finite  element  code,  UFPAC,  which 
was  developed  in  this  work.  Techniques  for  estimating  dynamic 
response  characteristics  of  the  structure  about  prestressed 

configurations  were  also  implemented  in  the  finite  element 
code . 
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The  finite  element  formulation  is  based  on  the 
incremental  virtual  work  principle  of  nonlinear  elasticity. 
The  equations  of  motion  are  linearized  about  the  prestressed 
configuration  so  that  the  nonlinearity  considered  in  the  final 
analysis  is  limited  to  the  stiffening/slackening  effect  of  a 
linear  preload.  The  required  response  characteristics  about 
the  prestressed  configuration  are  obtained  by  solving  the 
resulting  equations  of  motion  appropriately. 

Loss  factor,  a measure  of  damping  ability,  was  a key 
parameter  in  this  study.  Finite  element  method  based 
techniques  to  estimate  overall  loss  factor  based  on  the 
knowledge  of  the  loss  factors  of  the  constituent  materials 
that  form  the  composite  structure  are  discussed.  The 
modifications  required  in  the  damping  estimation  models  to 
account  for  the  prestressed  state  of  the  structure  are 
presented. 

A simple  experimental  method  to  measure  damping  ability 
is  described.  The  experiments  were  used  mainly  to  provide 
loss  factor  data  which  were  supplied  as  input  in  the  finite 
element  analysis.  Some  of  the  finite  element  predictions  were 
compared  with  experimental  results  to  verify  the  analytical 
model  for  the  prediction  of  damping.  Several  other  examples 
that  verify  the  accuracy  of  the  finite  element  formulation  are 
also  presented. 
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The  analysis  tools  developed  were  used  to  examine  the 
effect  of  viscoelastic  damping  treatments  on  the  response 
characteristics  of  laminated  composite  structures.  Parametric 
studies  which  show  the  effect  of  the  different  relevant 
factors  on  the  overall  damping  ability  are  presented.  The 
results  show  that  if  appropriately  designed,  it  is  possible  to 
construct  integrally  damped  composite  structural  elements  with 
much  improved  dynamic  response  characteristics. 
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CHAPTER  1 
INTRODUCTION 

1 . 1 Background 

Problems  associated  with  vibration  control  arise 
everywhere,  and  all  vibrating  systems  are  subject  to  some 
extent  of  damping.  Damping  is  the  dissipation  of  vibratory 
energy  from  a structure  executing  cyclic  motion.  In 
engineering  structures  which  encounter  dynamic  load 
environments,  it  is  almost  always  desirable  for  the  structure 
to  have  high  damping  ability.  Different  design  approaches  are 
used  for  structural  damping  improvement,  and  therefore, 
increase  in  the  durability  of  primary  and  secondary  structural 
components.  These  approaches  include  active  control,  passive 
control  and  remedial  control. 

Active  control  uses  external  devices  called  sensors  and 
actuators  to  shape  the  response  of  the  system.  In  general, 
active  control  uses  real  time  measurements  of  system  response. 
Passive  control,  on  the  other  hand,  uses  current  changes  in 
material  parameters  to  alter  damping  ability  and  system 
response.  Remedial  control  is  the  name  given  by  the  author  to 
situations  where  damping  ability  is  not  considered  as  a design 
parameter  and  is  added  to  the  system  as  dictated  by  the 
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problems  encountered  in-service.  In  advanced  aircraft  and 
aerospace  applications  damping  is  a key  design  parameter  and 
should  be  made  an  integral  part  of  system  design. 

With  the  use  of  appropriate  structural  and  damping 
materials,  it  is  possible  to  design  integrally  damped  members 
with  distributed  stiffness  and  damping  characteristics. 
However,  in  order  to  make  damping  a creative  force  in  design 
it  should  be  possible  to  create  and  control  the  magnitude  and 
type  of  damping.  Therefore,  analysis  techniques  are  necessary 
which  can  predict  the  damping  ability  of  the  structure 
accurately.  In  the  present  study,  the  concept  of  damping  by 
design  is  extended  only  to  tailoring  the  overall  passive 
damping  ability  of  the  structure. 

A number  of  damping  mechanisms  are  operative  in  a 
vibrating  structure  which  can  cause  energy  dissipation  and 
modify  its  response  characteristics.  A vibrating  structure 
always  interacts  with  the  surrounding  medium  and  may  dissipate 
energy  by  radiation  away  from  the  system.  For  example,  in 
acoustic  radiation  damping  vibratory  energy  is  radiated  away 
as  sound  [1].  Another  damping  mechanism  is  due  to  the 
frictional  forces  which  arise  when  contacting  surfaces  in  a 
structure  rub  against  each  other.  In  this  case,  the  vibratory 
energy  of  the  structure  is  converted  to  heat  energy.  The 
mechanisms  described  above  are  nonmaterial  damping  mechanisms. 
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After  all  such  mechanisms  are  accounted  for,  there  are 
still  a large  number  of  mechanisms  that  are  operative  within 
the  material  volume  which  can  cause  energy  dissipation  as  the 
material  deforms  [1-6].  The  different  damping  mechanisms  and 
the  conditions  under  which  they  are  active  are  listed  in 
Nashif , Jones  and  Henderson  [1],  All  these  mechanisms  are 
associated  with  internal  reconstructions  of  the  material  at 
different  scales  and  are  classified  under  material  or  passive 
damping.  Such  mechanisms  are  highly  nonlinear,  and  a general 
mathematical  description  is  quite  complicated  and  usually  not 
used  in  vibration  analysis.  Simpler  damping  models  which 
capture  the  effect  of  the  mechanism  on  system  response  have  to 
be  used. 

Of  the  different  damping  mechanisms,  viscoelastic 
damping,  which  is  exhibited  strongly  in  many  polymeric  and 
glassy  materials,  is  the  most  attractive  and  has  many 
possibilities  for  practical  applications  [7-10].  The  coiling 
and  uncoiling  of  the  three-dimensional  network  of  long  chain 
molecules  give  rise  to  the  high  damping  behavior  of  these 
materials  [11].  The  macromechanical  effect  of  this  molecular 
behavior  is  that  energy  dissipation  is  almost  entirely  due  to 
shear  dissipation  and  independent  of  dilatational  deformation 
[12].  It  is  possible  to  design  polymeric  materials  to  have 
high  damping  over  a wide  range  of  temperatures,  and  therefore, 
effective  damping  treatments  can  be  designed  from  the 
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Since  energy  dissipation  is  primarily  due  to  shearing  of 
the  damping  material,  the  damping  materials  should  be  included 
in  the  structure  in  a constrained  form.  The  damping  material 
may  be  included  as  an  integral  part  of  the  structure  or  may  be 
applied  to  the  surface  of  the  structure  as  dictated  by  the 
design  situation.  Surface  damping  treatments,  as  the  name 
suggests,  are  applied  to  the  surface  of  the  structure.  They 
are  classified  as  free-layer  damping  treatments  or 
constrained-layer  damping  treatments  depending  on  whether  the 
damping  layer  is  free  or  constrained  by  a stiff  constraining 
layer.  The  subject  of  viscoelastic  surface  damping  treatments 
and  the  mechanism  by  which  energy  is  dissipated  in  structures 
which  include  such  damping  treatments  are  discussed  in  greater 
detail  in  Chapter  2. 

Viscoelastic  materials  have  been  used  for  quite  some  time 
to  solve  a variety  of  noise  and  vibration  problems  associated 
with  sheet  metal  structures.  Consequently,  a large  amount  of 
work  has  been  done  in  the  area  of  analysis  of  damping 
treatments  applied  to  such  structures. 

Many  current  applications  require  the  use  of  anisotropic 
composite  materials  which  have  some  distinct  advantages  over 
conventional  monolithic  materials  [13-14].  The  improved 
specific  strength  and  stiffness  properties  of  composite 
materials  have  made  them  prime  candidates  in  several  advanced 
applications.  Composites  also  have  much  higher  inherent 
damping  ability  than  conventional  monolithic  materials,  and 
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their  damping  ability  can  be  tailored  to  a certain  degree.  In 
the  recent  past,  the  optimization  of  internal  material  damping 
of  various  types  of  advanced  composites  has  been  examined  [15- 
19].  These  analytical  and  experimental  investigations 
included  aligned-short-f iber  composites,  aligned-short-f iber 
offaxis  composites,  randomly  oriented  short-fiber  composites, 
and  two-  and  three-dimensional  analysis  of  laminated 
composites.  Parameters  such  as  loading  angle,  fiber  aspect 
ratio,  fiber  tip  spacing,  ply  thickness,  and  damping  ratio 
between  fiber  and  matrix  materials  were  adjusted  to  improve 
the  performance  of  composite  materials  in  a dynamic  load 
environment.  However,  an  increase  in  damping  ability  is 
usually  associated  with  a corresponding  stiffness  loss,  and 
this  limits  the  usable-feasible  options  to  increase  damping. 
Further  increase  in  damping  can  be  achieved  by  using 
viscoelastic  damping  materials  with  composite  structures,  just 
as  can  be  done  in  the  case  of  monolithic  structures. 

In  the  present  work  the  potential  for  using  viscoelastic 
materials  to  further  increase  the  damping  ability  of  laminated 
composite  structural  elements  is  examined.  Viscoelastic 
damping  treatments  and  the  development  of  a finite  element 
analysis  model  which  is  specially  suited  for  the  analysis  of 
viscoelastically  damped,  prestressed,  generally  laminated 
anisotropic  structures  are  described  in  the  remaining 
chapters . 
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1.2  Objective  and  Outline  of  This  Study 

The  analysis  of  vibrations  of  damped  structures  is 
extremely  important  in  designing  aircraft  and  aerospace 
structures.  As  was  discussed  in  the  previous  section,  one 
approach  to  improving  the  vibration  response  characteristics 
of  a structure  is  by  adding  viscoelastic  damping  layers.  In 
order  to  design  effective  passive  dampers,  accurate  analysis 
capability  is  reguired  to  predict  the  amount  of  damping.  In 
a typical  viscoelastically  damped  structure,  the  integral 
design  procedure  may  result  in  widely  dispersed  damping 
(viscoelastic  material  applied  only  to  some  portions  of  the 
structure)  , and  the  analysis  of  such  a structure  is  quite 
complicated. 

The  main  objective  of  this  research  was  to  develop  an 
efficient  finite  element  model  for  the  analysis  of 
viscoelastically  damped  composite  structural  elements. 
Emphasis  was  placed  on  the  efficiency  of  the  finite  element 
model,  since  three-dimensional  models  of  sandwich 
configurations  lead  to  excessively  large  degrees  of  freedom. 
Another  problem  of  interest  was  to  develop  the  ability  to 
model  generally  laminated  composites.  In  the  present  study, 
the  effect  of  prestress  on  the  vibration  damping 
characteristics  was  also  included.  The  consideration  of 
prestress  is  important  because  response  characteristics  may  be 
altered  significantly  by  prestress. 
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The  finite  elements  developed  were  implemented  in  a 
stand-alone  finite  element  analysis  program,  UFPAC  (Appendix 
C)  . The  finite  element  formulation  is  based  on  the  principle 
of  incremental  virtual  work  and  is  described  in  Chapter  3,  and 
the  element  matrix  formulations  are  described  in  Chapter  4. 
The  techniques  for  estimation  of  the  dynamic  response 
characteristics  of  a structure  from  the  finite  element 
equations  of  motion  are  described  in  Chapter  5.  Finally,  the 
analysis  tools  are  used  to  examine  the  improvement  in  the 
damping  performance  of  composite  structural  elements  that  can 
be  achieved  by  using  viscoelastic  dampers. 


CHAPTER  2 

VISCOELASTIC  DAMPING  TREATMENTS 

Viscoelastic  treatments  can  be  used  to  solve  a variety  of 
resonant  noise  and  vibration  problems.  Such  treatments  may  be 
integrated  into  a structural  component  or  may  be  applied  to 
the  surface  of  an  existing  structure  to  provide  high  damping 
ability  over  a wide  temperature  and  freguency  range.  Surface 
damping  treatments  are  commonly  classified  under  one  of  two 
categories,  extensional  or  shear,  based  on  the  deformation 
experienced  by  the  viscoelastic  damping  layer.  Extensional 
damping  treatments  are  also  called  free-layer  damping 
treatments  simply  because  one  surface  of  the  damping  layer  is 
free  or  unconstrained.  In  this  type  of  damping  treatment, 
energy  is  dissipated  as  the  vibratory  energy  is  transferred 
from  the  structure  to  the  damping  layer  which  experiences  the 
surface  deformation  of  the  base  structure.  As  should  be 
expected,  damping  ability  increases  with  increase  in  the 
thickness  of  the  damping  layer  [1,20].  However,  use  of  a 
thick  layer  of  damping  material  has  an  associated  weight 
penalty,  and  the  redistributed  mass  and  stiffness 
characteristics  may  alter  response  characteristics  more  than 
is  permissible. 

Viscoelastic  damping  treatments  are  most  effective  if  the 
damping  material  is  included  in  a constrained  form.  By 
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constraining  the  viscoelastic  layer  so  that  a large  difference 
in  stiffness/rigidity  is  present  across  the  viscoelastic 
damping  layer,  it  is  possible  to  induce  shear  deformation  in 
the  damping  layer.  Such  a design  will  yield  the  best  results 
for  a given  volume  of  the  damping  material,  since  most  real 
materials  are  found  to  be  viscoelastic  in  shear  and  elastic  in 
dilatation  [1].  In  the  case  of  surface  shear  damping 
treatments,  the  damping  layer  is  forced  to  deform  in  shear  by 
constraining  the  soft  viscoelastic  layer  by  a thin,  stiff 
constraining  layer.  As  the  structure  deforms  cyclically,  the 
stiff  constraining  layer  will  force  the  damping  layer  to 
deform  in  shear.  The  concept  is  illustrated  in  Figure  2.1. 

All  real  elastomeric  materials  have  very  strong 
temperature  and  freguency  dependent  material  properties.  The 
effectiveness  of  a damping  treatment  depends  on  the  conditions 
under  which  it  has  to  operate.  The  effect  of  temperature  on 
the  properties  of  a typical  viscoelastic  material  is  shown  in 
Figure  2.2.  At  low  temperatures,  the  stiffness  of  the 
material  is  high  and  the  damping  ability  is  low.  In  this 
situation,  the  energy  dissipation  ability  of  the  material  is 
low  because  of  the  low  damping  ability  and  the  reduced  shear 
deformation  of  the  viscoelastic  material.  Reduced  shear 
deformation  is  caused  by  the  rigid  coupling  between  the 
viscoelastic  layer  and  the  constraining  layer.  At  higher 
temperatures,  where  the  material  is  in  the  soft  rubbery 
region,  the  dissipation  is  again  small  because  the  loss 
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Figure  2.1.  Viscoelastic  Surface  Damping  Treatment 
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Figure  2.2.  Variation  of  Material  Properties  with  Temperature 


12 


modulus  of  the  damping  material  is  low.  Somewhere  in  between, 
in  the  glass  transition  region,  the  material  possesses  an 
optimal  loss  modulus  so  that  the  damping  ability  is  a maximum. 
It  is  always  desirable  to  have  operating  conditions  in  this 
zone.  Several  damping  materials  are  available  which  have 
transition  regions  over  a range  of  temperatures  from  -200  F 
to  2000  F [21-22].  A large  amount  of  work  has  been  done  in 
the  area  of  characterization  of  material  properties  of  such 
materials.  The  properties  as  a function  of  frequency  and 
temperature  are  presented  concisely  using  a reduced  frequency 
nomogram  [23].  The  reduced  frequency  chart  for  the  damping 
tape  ISD1121  is  given  in  Figure  2.3.  This  damping  material 
was  used  in  all  the  problems  discussed  in  the  results  section. 

Fundamental  work  in  the  area  of  the  analysis  of  sandwich 
viscoelastic  beams  was  done  by  Kerwin  [24],  Kerwin  considered 
a beam  coated  by  an  infinitely  long  constrained  damping 
treatment.  Kerwin' s work  establishes  a complex  stiffness  for 
the  beam  assuming  that  (a)  a perfect  bond  exists  between  the 
base  structure  and  the  damping  layer,  (b)  a major  part  of  the 
damping  is  due  to  the  shearing  of  the  viscoelastic  layer,  (c) 
the  deflection  of  the  beam  does  not  vary  through  the 
thickness,  and  (d)  the  beam  is  simply  supported  or 
infinitely  long  so  that  beam  fixity  may  be  neglected.  Kerwin 
showed  that  the  normal  force  in  the  viscoelastic  layer  is 

1A  class  of  room  temperature  damping  tape  manufactured  by 
3M  company. 
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small  in  comparison  to  the  force  in  the  base  beam  and  the 
constraining  layer  and  can  be  neglected.  This  assumption 
reduces  the  complexity  of  analysis  significantly  and  has  been 
used  repeatedly  in  the  analysis  of  conventional  sandwich  beams 
and  beams  treated  by  a constrained  viscoelastic  damping 
layer.  Even  though  Kerwin's  analysis  is  valid  only  for 
simply  supported  beams,  it  has  been  used  for  beams  with  other 
boundary  conditions  by  patching  on  a correction  factor  [25]. 

Ditaranto  [26]  developed  a more  general  6th  order  theory 
for  beams  with  arbitrary  boundary  conditions  by  using  a more 
descriptive  set  of  coordinates.  Ditaranto  presupposed  a 
solution  of  a certain  form  which  led  to  the  deduction  that  the 
relation  between  frequency  and  loss  factor  was  independent  of 
boundary  condition.  Mead  and  Markus  [27-28]  reexamined  the 
assumption  of  Ditaranto  and  concluded  that  it  was  without 
basis.  They  also  showed  that  maximum  loss  factors  of  beams 
with  different  boundary  conditions  occur  at  different 
frequencies.  Exact  fundamental  frequencies,  loss  factors  and 
optimal  design  of  beams  with  different  boundary  conditions 
have  also  been  reported  in  other  works  [29,30].  In  each  case 
the  analysis  presented  is  suitable  for  conventional  sandwich 
beams  and  beams  fully  treated  by  constrained  damping 
treatments . 

Parfitt  [31]  examined  the  effect  of  cuts  on  the  damping 
ability  of  a constrained  layer  damping  treatment.  He  showed 
that  the  damping  ability  can  be  increased  by  periodic 
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interruption  or  segmenting  of  the  constraining  layer. 
Segmenting  induces  shear  in  the  viscoelastic  damping  layer  and 
thereby  increases  damping  ability.  Lazan  [32]  showed  that  the 
amount  of  damping  can  be  increased  by  using  an  alternately 
anchored  multiple  layer  treatment. 

Plunkett  [33]  showed  that  when  properly  assembled,  the 
multiple  layers  need  not  be  anchored.  Plunkett  also  presented 
a theory  for  finite  length  and  thickness  of  treatment  and 
compared  the  predictions  of  the  theory  with  experimental 
results.  Additional  simplifying  assumptions  were  made  by 
Plunkett  since  his  theory  was  specialized  for  the  analysis  of 
beams  treated  by  constrained  damping  treatments.  Plunkett's 
assumptions  are  characteristic  of  similar  efforts  and  are 
listed  below.  (a)  The  thicknesses  of  the  constraining  layer 
and  the  damping  layer  are  very  small  when  compared  to  that  of 
the  base  structure;  thus  the  bending  effects  in  these  layers 
are  negligible.  This  implies  that  the  constraining  layer  is 
subjected  to  tension  only,  and  the  constrained  layer  is 
subjected  to  shear  only  (Kerwin's  assumption).  (b)  The 
damping  layer  is  linearly  viscoelastic,  (c)  The  constraining 
layer  is  purely  elastic,  (d)  The  shear  strain  is  uniform 
through  the  thickness  of  the  viscoelastic  layer,  (e)  The 
normal  stress  is  uniform  through  the  thickness  of  the 
constraining  layer.  (f)  The  elastic  moduli  of  the 
viscoelastic  layer  are  small  in  comparison  with  those  of  the 
base  structure  and  constraining  layer. 
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The  analysis  procedures  discussed  above  are  useful  for 
examining  the  basic  mechanisms  of  surface  damping  treatments 
and  are  applicable  to  problems  involving  specialized  geometry, 
loadings  and  boundary  conditions.  Even  for  simple  geometries, 
algebraic  or  numerical  solution  of  the  problem  tends  to  be 
complicated  and  lengthy.  In  any  case,  these  methods  are  not 
usable  when  the  base  beam  or  plate  is  a generally  laminated 
composite.  Practical  problems  usually  involve  complicated 
structural  geometries  and  boundary  conditions,  and  only  a 
portion  of  the  structure  has  damping  treatment  applied  to  it. 
Since  much  of  the  difficulty  of  analyzing  constrained  layer 
damping  treatment  is  due  to  complicated  geometries,  it  is 
normal,  as  in  the  case  of  undamped  structures,  to  look  at 
finite  element  technigues  for  solutions  to  the  problem.  By 
using  the  finite  element  method,  arbitrary  boundary  conditions 
and  loadings  can  be  modeled  guite  easily.  The  finite  element 
method  is  currently  popular  for  analyzing  and  designing  damped 
structures . 

The  works  related  to  finite  element  analysis  of 
viscoelastically  damped  structures  deal  mainly  with  approaches 
to  damped  structural  design  in  the  context  of  implementation 
using  existing  general  purpose  finite  element  analysis  codes 
[34-38].  Soni  and  Bogner  [39,40]  modified  a general  purpose 
finite  element  analysis  code,  MAGNA,  for  the  analysis  of 
damped  structures.  MAGNA  is  strongly  oriented  toward 
isotropic  elastic  solid  continua.  The  analytical  methods  to 
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predict  damping  of  constrained-layer  damped  structures  have 
been  evaluated  in  several  papers  [34-38,41],  and  are  discussed 
in  Chapter  5 . 

In  the  present  work,  structural  finite  elements  which  are 
specially  suited  to  model  generally  laminated,  prestressed 
composite  structures  subjected  to  constrained  viscoelastic 
damping  are  developed.  An  efficient  modification  of  the  modal 
strain  energy  approach  [42]  for  the  estimation  of  damping 
ability  about  prestressed  configurations  is  also  presented. 
The  finite  elements  and  the  damping  estimation  procedures  are 
implemented  in  a program,  UFPAC,  developed  by  the  author. 
Details  about  the  program  and  its  execution  are  presented  in 
Appendix  C. 


CHAPTER  3 

FINITE  ELEMENT  FORMULATION 


In  this  chapter,  the  finite  element  equations  for  the 
analysis  of  a structure  about  a linearly  prestressed 
configuration  are  developed  based  on  the  principle  of 
incremental  virtual  work  [43-49].  The  incremental  virtual 
work  principle  is  expressed  in  terms  of  incremental  kinematic 
and  force  quantities,  and  is  derived  formally  from  the 
principle  of  virtual  work.  The  principle  of  virtual  work, 
also  called  the  principle  of  virtual  displacements,  is  an 
equivalent  approach  to  establishing  the  equilibrium 
expressions  of  a structure.  It  forms  the  basis  of  most 
existing  general  purpose  finite  element  codes  due  to  the 
relative  simplicity  with  which  the  formulation  can  be 
implemented. 

According  to  this  principle  [50],  if  Jv  T^  Se13  dV  is 
equal  to  the  virtual  work  of  all  external  forces  for  all 
kinematically  admissible  virtual  displacement  fields 
(compatible  virtual  displacements  satisfying  all  essential 
boundary  conditions),  then  the  stress  field,  Tid , satisfies  all 
equilibrium  requirements.  The  principle  of  virtual  work  is 
written  in  indicial  notation  as 
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f{T2jdeij  + P dv  = J fibui  dv  + f t i6ui  dA 

v v0  s 


T-  • 
ij 


<5e. 


ij 


fi 

ti 

Ui 


Cauchy  stress  tensor 

infinitesimal  strain  tensor  calculated  from  <5^ 
body  forces 
surface  tractions 

infinitesimal  displacements  from  equilibrium 
configurations 


accelerations 

Forces  due  to  the  state  of  motion  of  the  body  are  accounted 
for  in  the  D'Alembert  sense.  The  finite  element  equilibrium 
equations  for  a displacement  based  linear  structural  analysis 
can  be  derived  by  appropriate  discretization  of  the  above 
virtual  work  expression.  The  kinematic  and  force  quantities 
are  defined  in  the  deformed  configuration,  but  the  distinction 
between  reference  states  is  usually  not  made  since  the 
displacements  and  strains  are  assumed  infinitesimal. 

Even  under  the  assumption  of  small  motion,  there  are 
certain  loading  situations  which  cannot  be  modeled  by  using 
the  linear  formulation.  One  such  situation  is  when  in-plane 
and  out-of-plane  loads  are  applied  simultaneously.  The  in- 


plane loads  will  affect  the  transverse  load  carrying  ability 
of  a structural  element,  and  the  coupling  between  inplane  and 
out  of  plane  effects  is  a nonlinear  effect.  The  principle  of 
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superposition,  which  holds  in  all  linear  analyses  will  not 
apply  for  such  loading  conditions  since  the  displacements  are 
not  a linear  function  of  the  load.  In  future,  the  term 
"prestress  effect"  will  be  used  to  define  the  coupling  between 
in-plane  and  out-of-plane  loads. 

The  continuum  equations  which  include  the  prestress 
effect  are  developed  from  the  incremental  virtual  work 
principle  of  nonlinear  elasticity  [43-49].  These  equations 
are  valid  for  a general  nonlinear  problem  and  include  all 
kinematic  and  geometric  nonlinear  effects.  The  complete 
formulation  is  presented  to  discuss  the  different  nonlinear 
effects;  however,  in  the  finite  element  code,  UFPAC,  that  was 
developed  in  this  work,  only  the  prestress  effect  was 
considered. 

The  totally  Lagrangian  formulation  for  the  analysis  of 
nonlinear  motion  is  discussed;  however,  other  formulations  are 
also  applicable  to  solution  of  such  problems  [43,46].  The 
choice  of  formulation  was  dictated  by  the  category  of  analysis 
that  was  of  interest  in  this  research.  In  a totally 
Lagrangian  formulation,  as  indicated  by  the  name,  all 
kinematic  and  force  quantities  are  referred  to  the 
referential/undeformed  coordinate  system.  The  virtual  work 
expression,  therefore,  has  to  be  written  in  terms  of 
quantities  that  are  defined  in  the  undeformed  configuration. 
Derivation  of  the  virtual  work  principle  in  the  reference 
configuration  can  be  found  in  Malvern's  work  [50]. 
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In  a problem  involving  finite  deformations,  the  final 
configuration  of  the  structure  is  unknown.  This  is  in 
contrast  to  infinitesimal  deformation  problems  where  it  can  be 
assumed  that  the  deformed  body  is  almost  similar  in 
configuration  to  the  undeformed  body.  In  the  incremental 
approach,  a seguence  of  motion  is  considered,  and  the 
incremental  displacement  and  force  quantities  are  calculated 
so  that  the  equilibrium  relations  are  satisfied  in  the  current 
configuration.  Tangent  measures  of  the  stiffness  matrix  and 
the  mass  matrix  are  evaluated  at  each  step.  Since  the 
configuration  of  the  structure  in  finite  deformation  problems 
can  change  considerably,  a stiffness  measure  based  on  initial 
undeformed  coordinates  will  not  be  valid  for  the  entire 
deformation  history.  This  is  the  primary  reason  why  a linear 
analysis  has  to  be  modified  in  geometrically  nonlinear 
problems . 

The  ease  with  which  a nonlinear  analysis  can  be 
implemented  depends  upon  the  stress  and  strain  measures  used. 
The  Cauchy  stress  tensor  is  defined  in  the  unknown  deformed 
configuration.  Also,  the  Cauchy  stress  tensor  is  not  rotation 
invariant.  The  Cauchy  stress  in  the  current  configuration,  in 
general,  cannot  be  obtained  by  simply  summing  the  incremental 
stress  components  due  to  increments  in  strains,  since  rigid 
body  rotation  will  change  the  components  of  the  Cauchy  stress 
tensor.  The  second  (symmetric)  Piola-Kirchhof f stress  tensor 
is  essentially  rotation  invariant  for  the  small  strain,  large 
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deformation  case  and  can  be  used  with  appropriate 
constitutive  relations  to  develop  an  efficient  incremental 
formulation. 

The  motion  of  a generic  three-dimensional  body  is  shown 
in  Figure  3.1.  The  motion  is  not  necessarily  time  dependent, 
but  the  time  variable  is  used  to  describe  the  loading  and 
motion  history  of  the  body.  The  objective  is  to  calculate  the 
final  equilibrium  configuration  of  the  body  at  the  end  of  the 
load  history  curve.  At  each  stage,  it  is  assumed  that  the 
solution  to  the  problem  at  discrete  time  points  from  time  t=0 
to  t=t  is  known  and  that  the  solution  at  time  (t+At)  is  to  be 
calculated. 

Consider  the  body  of  Figure  3.1,  whose  initial 
configuration  is  denoted  by  C0  and  in  which  the  Cartesian 
coordinates  are  assigned  to  a point  in  the  structure  in  C0. 
After  subsequent  deformation  of  the  body,  the  position  of  the 
same  particle  is  given  by  xL  in  its  current  configuration  C1. 
This  state  is  an  intermediate  state  in  the  typical  sequence  of 
motion  described  above.  State  C2  is  the  deformed  state,  the 
equilibrium  configuration  at  time  (t+At) , to  be  determined 
after  the  current  increment  of  load  is  applied.  The 
configuration  in  the  deformed  state  is  actually  calculated  by 
evaluating  the  incremental  solution  between  the  states  C2  and 
C2  and  updating  the  intermediate,  Cx  state  deformation.  In 
other  words,  the  solution  procedure  involves  the  summing  of 
displacement  increments  due  to  corresponding  load  steps  to  the 
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Intermediate 

State 


Figure  3.1.  Successive  Deformation  States 
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point  where  the  out-of-balance  load  is  zero.  The  out-of- 
balance  load  is  the  portion  of  the  total  load  not  balanced  by 
internal  stresses  that  are  developed  in  the  continuously 
deforming  structure. 

The  principle  of  virtual  work  [50]  in  state  a,  in  terms 
of  tensors  referred  to  the  undeformed  state,  is  written  using 
indicial  notation  as 

/{  AAA  + oP  «0AUi}  dV  = 

(3.2) 

/ afAui  dV  + f atd8 aud  dS 

vo  dV0 

where 


sij 

second  symmetric  Piola-Kirchhof f stress  tensor 

Eij 

Green-Lagrange  strain  tensor 

fi 

body  forces 

ti 

surface  tractions 

Ui 

displacements 

Ui 

accelerations 

The  left 
quantity 
variables 
The  force 
forces  in 

subscript  denotes  the  configuration  in  which  a 
is  being  defined,  and  all  kinematic  and  force 
are  referred  to  the  initial  undeformed  geometry, 
due  to  the  state  of  motion  is  included  as  body 
the  D'Alembert  sense. 

Displacements , ^Ui , are  resolved  in  the  directions  of  Xx 
and  the  strains  are  given  by  the  Green-Lagrange  strain  tensor. 
The  Green-Lagrange  strain  tensor  in  state  a is  given  by 
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(3.3) 


The  stresses  are  the  Cartesian  components  of  the  second 
(symmetric)  Piola-Kirchhof f stress  tensor  [50].  The  second 
Piola-Kirchhof f stress  gives  a fictitious  force  in  the 
reference  configuration,  related  to  the  force  in  the  current 
configuration  in  the  same  way  as  a material  line  in  the 
referential  configuration  is  related  to  the  corresponding  line 
in  the  current  configuration.  The  components  of  the  second 
Piola-Kirchhof f in  state  a are  defined  in  terms  of  Cauchy 
stress  tensor  in  state  a as 
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or  in  matrix  form 
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where  F is  the  deformation  gradient  tensor,  whose  rectangular 
Cartesian  components  are 
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The  ratio  of  the  densities  can  be  expressed  in  terms  of  the 
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determinant  of  the  deformation  gradient  matrix.  If  the  Cauchy 
stress  tensor  and  the  deformation  gradient  are  known,  the 
Piola-Kirchhof f stress  tensor  can  be  calculated  easily. 
Similarly,  the  Cauchy  stress  tensor  can  also  be  calculated 
from  the  Piola-Kirchhof f stress  tensor.  Since  the  Cauchy 
stress  tensor  is  symmetric,  it  follows  that  the  Piola- 
Kirchhof  f stress  tensor  is  also  symmetric.  The  important 
implication  of  this  fact  is  that  the  finite  element  stiffness 
matrix  of  the  structure  will  always  be  symmetric. 

The  incremental  deformation  between  states  Cx  and  C2  is 
given  by  displacement  uir  where 

ui  = 2 ui  ~ iui  (3.7) 


All  quantities  with  no  left  subscripts  are  incremental 
quantities.  The  incremental  strain  is  derived  by  subtracting 
the  expression  for  strain  at  state  1 from  that  at  state  2,  and 
is  given  by 


ij  2^ij  1^. 


ij 


2 (uti  + Uj.i  + 1 Uk,i  luk.j  + luk.juk.i) 


+ 


Uk  .j) 


(3.8) 


= eij  + ci7- 

where  e^  and  C ij  are  the  linear  and  the  nonlinear  parts  of  the 
incremental  strain  respectively.  All  quantities  with  no  left 
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subscripts  represent  incremental  variables.  Incremental 
stress  is  represented  similarly  as  the  difference  between  the 
stress  at  states  1 and  2 respectively. 

Eij  ~ 2 Eij  ~ i Eij  (3.9) 

The  statement  of  incremental  virtual  work  is  derived  by 
subtracting  the  virtual  work  expressions  for  state  1 from  the 
expression  for  state  2.  The  expression  for  virtual  work  in 
state  2 in  terms  of  incremental  strains  and  stresses,  and  the 
known  stresses  and  strains  at  state  1 is 


/{(  i sij  + sij)  Mij  + oP(  i Ui  + Ui)*Ui}  dV  = 

V0 

(3 . 10) 

/(  ifi  + dV  + [(  1ti  + dS 

vo  dv0 


By  subtracting  Eguation  3.2  from  Equation  3.10,  and 
recognizing  that  62Eij=5Eij,  the  incremental  virtual  work 
principle,  in  terms  of  the  incremental  variables  is 


[SijbEij  dV  + f 1Sij  5£ij  dV  + f 0pui  dV  = 

Vo  v0 

(3.11) 

f f^UidV+J  ti6ui  dS  - f xSi:j  be±j  dV 
v0  dv0  V0 

By  assuming  a linear  relation  between  the  increments  of  Piola- 
Kirchhoff  stress  and  Green-Lagrange  strain,  S1J=CijrsErs , the 
statement  of  incremental  virtual  work  can  be  expressed 
entirely  in  terms  of  displacement  variables  as 
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f Cljr sE'JEu  dV  + f ±Si:j  b^  + / oPii,  dV  = 
v0  V0  v0 

f fi  8ui  dV'  + f ti6ui  dS  - f 6e ±j  dV 

v0  9V0  V0 


To  obtain  an  efficient  solution  procedure  to  the  nonlinear 
equation,  it  is  beneficial  to  linearize  the  incremental 
virtual  work  expression.  The  equation  is  linearized  by 
assuminq  linear  expressions  for  nonlinear  incremental  strains. 
The  linearized  virtual  work  equality  is 


/ Cijrsezsbei:j  dV  + f ^ SCij  dV  0pui  dV  = 

v0  V0 
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Equation  3.13  is  simplified  further  to  include  just  the 
nonlinear  effects  due  to  a linear  preload.  This  effect  has 
been  identified  earlier  as  the  prestress  effect.  In  all 
problems  considered  in  the  present  development,  it  is  assumed 
that  the  final  configuration  attained  by  the  structure  is 
reached  throuqh  a sequence  of  two  different  motions.  The 
first  increment  is  due  to  motion  induced  by  a known  initial 
load  followed  by  the  second,  a steady  state  harmonic  motion 
superposed  upon  this  initial  deformed  state.  Another  problem 
that  is  of  interest  is  the  calculation  of  the  modal  parameters 
about  the  prestressed  configuration.  The  first  increment  is 
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assumed  linear,  and  the  stiffness  characteristics  of  the 
structure  are  linearized  about  the  prestressed  configuration. 

The  sequence  of  motions  of  the  generic  body  shown  in 
Figure  3.1  can  be  used  to  define  the  problem  as  follows. 
Consider  a body  whose  initial  configuration  is  denoted  by  C0 
and  in  which  Cartesian  coordinates,  XA,  are  assigned  to  a 
point  in  the  structure  in  the  undeformed  configuration.  After 
subsequent  deformation  of  the  body,  due  to  a known  static 
linear  preload,  the  same  particle  occupies  the  location  xi  in 
the  intermediate  configuration  C1.  C2  is  the  final 
configuration  to  be  calculated,  due  to  second  increment  of  the 
load.  Since  only  small  deformations  are  considered  in  both 
motions,  a further  assumption  is  made  that  the  volume  of  the 
body  does  not  change  so  that  the  Piola-Kirchhof  f stress 
measure  and  the  Cauchy  stress  measure  are  approximately  the 
same.  Therefore,  the  Cauchy  stress  at  the  prestressed  state 
is  used  to  calculate  the  contribution  of  the  nonlinear 
prestress  effect  to  the  overall  stiffness  matrix.  The 
required  response  characteristics  of  the  structure  can  be 
calculated  by  using  the  techniques  described  in  Chapter  5. 

A family  of  laminated  structural  and  continuum  elements 
is  developed  from  the  general  formulation  presented,  specially 
suited  for  the  analysis  of  viscoelastically  damped  composite 
structures.  Only  structures  treated  with  constrained 
viscoelastic  surface  damping  treatments  are  considered.  As 
described  in  a previous  section,  most  of  the  energy 
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dissipation  in  such  structures  is  due  to  the  shearing  of  the 
viscoelastic  core.  Therefore,  it  is  required  that  the  finite 
element  model  used  should  be  able  to  accurately  model  the  core 
shearing.  Special  emphasis  is  also  placed  on  the  efficiency 
of  the  model,  since  a fully  three  dimensional  model  for  a 
laminated  sandwich  can  easily  run  in  to  large  problem  sizes. 
In  subsequent  chapters,  the  formulation  of  specialized  plate 
and  beam  elements,  and  the  finite  element  analysis  of 
laminated  plates  and  beams  treated  with  constrained 
viscoelastic  damping  layers  are  presented. 


CHAPTER  4 

DERIVATION  OF  ELEMENT  MATRICES 


The  finite  elements  formulated  in  this  section  are 
specially  suited  to  model  three  layer  sandwich  structures. 
The  term  "sandwich  structure"  is  usually  used  to  define  a 
structure  with  a low  weight  core  and  high  stiffness  facings. 
In  the  present  context,  the  term  is  used  to  define  a stiff 
base  treated  by  a constrained  damping  treatment  (Figure  2.1). 
In  such  an  assembly  the  load  bearing  ability  is  retained 
primarily  in  the  base  structure,  and  the  purpose  of  the 
applied  damping  material  is  purely  to  increase  the  passive 
damping  ability  of  the  structure. 

The  base  structure  and  the  constraining  layer  are  modeled 
with  beam  or  plate  elements,  and  the  viscoelastic  core  is 
modeled  with  two-dimensional  plane,  or  three-dimensional  solid 
elements.  One  feature  of  the  formulation  is  that  the  plate 
and  beam  elements  are  naturally  compatible  with  the  continuum 
elements.  Plate  elements  which  are  formulated  conventionally 
(with  reference  plane  as  mid-plane)  are  incompatible  with 
continuum  elements,  and  special  constraints  have  to  be  used  to 
ensure  compatibility. 

In  the  following  sections,  the  derivation  of  the  tangent 
stiffness  matrix,  the  geometric  stiffness  matrix  and  the  mass 
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matrix  of  the  structural  and  continuum  elements  is  presented. 
The  plate  and  solid  elements  are  isoparametric  elements,  a 
concept  attributed  to  Taig  [51]  and  Irons  [52].  A key  feature 
of  the  isoparametric  formulation  is  that  the  same 
interpolation  functions  are  used  to  define  the  coordinates  and 
the  displacements. 

4.1  Laminated  Variable-Node  Offset  Plate  Element 
The  element  can  be  described  as  a variable-eight-node, 
shear-deformable,  laminated,  offset  plate  element.  The  plate 
element  in  the  element  coordinate  system,  and  the  active 
degrees  of  freedom  per  node  are  shown  in  Figure  4.1.  The 
nodes  numbered  5 through  8 are  variable,  and  may  be  added  to 
permit  quadratic  behavior  along  an  edge.  When  just  nodes  1 
through  4 are  present  the  element  behaves  like  a bilinear 
quadrilateral  plate,  and  when  all  eight  nodes  are  present  the 
element  exhibits  quadratic  behavior.  This  flexibility  allows 
for  simple  transitions  from  a coarse  to  a fine  mesh,  and  also 
the  use  of  higher  order  interpolation  functions  where 
required.  The  term  shear  deformable  implies  that  the 
kinematic  assumptions  are  based  on  the  Reissner-Mindlin  plate 
theory.  The  fiber  rotation  and  deflection  are  kinematically 
independent,  and  this  is  a departure  from  the  classical 
Poisson-Kirchhof f plate  theory.  The  element  may  either  be 
monolithic  or  laminated,  with  properties  varying  through  the 
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Sign  convention  for  rotations 


02 


Figure  4.1.  Laminated  Offset  Plate  Element 
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thickness.  Finally,  in  conventional  plate  elements,  the 
reference  plane  in  which  the  displacements  and  rotations  are 
defined  is  also  the  geometric  midplane.  In  the  present 
formulation,  the  reference  plane  may  be  located  at  any  plane 
along  the  thickness  dimension.  This  can  be  done  by  including 
an  offset  parameter  in  the  formulation  which  defines  the 
perpendicular  distance  between  the  midplane  and  the  reference 
plane . 

The  elements  based  on  the  Reissner-Mindlin  theory  require 
only  C°  continuity,  and  permits  the  use  of  a variety  of 
interpolation  schemes.  Also,  the  problems  associated  with 
elements  based  on  the  classical  theory,  which  requires  C1 
continuity,  are  avoided  [53]. 

4.1.1  Elastic  Stiffness  Matrix 

The  main  assumptions  of  the  Reissner-Mindlin  plate 
theory  are  (a)  azz= 0,  (b)  plane  sections  remain  plane  and,  (c) 
deflection  does  not  vary  through  the  thickness.  The  assumed 
displacement  field  is 

u(x,y,z)  = u0(x,y)  - zO^x^) 

v(x,y,z)  = v0(x,y)  - z02(x,y)  (4.1) 

w(x,y,z)  = w (x, y) 

where 

u displacement  in  the  positive  x direction 

v displacement  in  the  positive  y direction 
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w displacement  in  the  positive  z direction 

0l7  02  rotations  (Figure  4.1) 

u0,  v0  displacements  in  the  plane  of  definition  of 
the  nodes  (reference  plane) 

The  stress-strain  relation  for  each  ply  with  respect  to  the 
orthotropic  axis  (axes  oriented  along  and  perpendicular  to  the 
fiber  axes  in  fiber  reinforced  composites)  is  [54] 
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The  stress-strain  relation  with  respect  to  the  laminate  axes 
oriented  at  an  angle  with  respect  to  the  ply  orthotropic  axes 
in  the  plane  of  the  laminate  is 
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The  eguation  is  rearranged  as 


[0]  [o]  t; 

[0]  [Qs]  N 


(4.3) 


(4.4) 


[Q]  and  [Qs]  are  the  transformed  ply  stiffnesses,  and  a1  and 
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cts  stand  for  inplane  stress  and  transverse  shear  stress 
respectively.  From  the  assumed  displacement  field.  Equation 
4.1,  the  strain-displacement  relations  are 


V 

^0,X 

0l,x 

ey 

V0,y 

02,v 

► = < 

U0.y+V0.X 

• - Z 1 

01,y+02,x 

vO,x-0l 

0 

eyz. 

W0.y~Q2 

0 

The  tangent  stiffness  matrix  for  the  plate  element  is 
calculated  from  the  first  term  in  the  linearized  expression  of 
incremental  virtual  work.  By  using  Equations  4.4  and  4.5,  the 
first  term  from  Equation  3.13  can  be  written  in  matrix  form  as 


f cijzsers^>eij  dv  = J{E}T[G]  {E\  dA 
v0  A 


where  the  vector  {E}  is 


{£} 


u, 


o,x 


V0.y 

( U0,y+V0,X ) 
(^o  ,y-02) 

01,  X 

02,  y 

(0l,y+02,x> 


(4.6) 


(4.7) 


and  the  matrix  [G]  is 
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[G] 


[A]  [ B ] [0] 

[B]  [ D ] [0] 

[0]  [0]  [F] 


(4.8) 


The  submatrices  [A],  [B]  , [D]  and  [F],  of  the  combined 

stiffness  matrix,  are  as  defined  below. 


z0+h/2 

Aij  = f [Q]  dz  i,  j = 1,2,3 

za-h/2 


z0*h/2 

Bij  = ~ f z [Q  ] dz  i,j  = 1,2,3 

z0-h/2 

(4.9) 

z0+h/2 

Dij  = f z2[Q  ] dz  i,j  =1,2,3 

z0-h/2 

z0*h/2 

FU  = f ti?s]  dz  i/j  =1,2 

z0-h/2 


z0  is  the  distance  between  the  midplane  of  the  laminate  and 
the  x-y  reference  plane,  and  h is  the  thickness  of  the 
laminate.  It  should  be  recognized  that  the  [A],  [B]  and  [D] 
matrices  are  equivalent  to  the  submatrices  of  the  combined 
stiffness  matrix  encountered  in  classical  laminate  analysis 
[55,56].  The  only  difference  is  that  the  stiffness 
coefficients  are  calculated  for  an  arbitrary  location  of  the 
reference  plane  in  the  present  discussion. 
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Within  each  element,  the  displacement  field  is  assumed  to 
be  of  the  form 

{u}  = [ N ] {U}e  (4.10) 


where  [N]  is  the  matrix  of  interpolation  functions,  <pir  and 
{U}e  is  the  vector  of  elemental  degrees  of  freedom.  The 
nonzero  entries  of  [N]  are 


Nu  = 4>* 


(4.11) 


(k=l, 6,  (j=  1,8  (i=( j-1) *6  + 1 ))) 

The  interpolation  functions  for  the  variable  eight-node 
master  plate  element  are 


<K 

({>2 

*3 

<t>4 


|(1+T1)  (1+0  - |4>5 


~ ( 1 t) ) (i+O  - -|<t>5 

|(1-T!)  (1-0  - |<t>6 

Id+t!)  (1-0  - |d)7 


4>5  = j (1-n2)  (i+O 
4>6  = ^U-r))  (1-0) 


ih 

4*7 


(4.12) 


4>7  = I (1-n2)  (i-O 


(j>8  = 4 d+n ) (1-0 


0i  = 0 if  node  i is  not  present  for  i=5,6,7,8 
By  using  the  expression  for  the  assumed  displacement  field, 
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Equation  4.10,  the  vector  {E}  of  Equation  4.7  can  be  expressed 
as 


{£}  = [P]  [N]T  {U}e  = [Q]{U)e  (4.13) 

where  [P]  is  a matrix  of  partial  differential  operators.  The 
differentiation  is  performed  with  respect  to  the  global 
coordinates.  The  operators  in  [P]  are 


[P] 
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dx 
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dy 
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dy 
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dx 
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ay 


_d_ 

dy 

_a_ 

dx 


(4.14) 


The  matrix  [Q]  is  the  matrix  of  derivatives  of  interpolation 
functions  formed  by  the  differential  operators  in  [P] 
operating  on  the  interpolation  functions  in  [N]  . By  using  the 
relation  in  Equation  4.13,  Equation  4.6  can  be  written  as 
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/ Cijrsers^>eij  dV  = J{6£)r[G]{£j  dA 
V0  A 

= |{6D}eT[^]  T[P]  r[G]  [P]  [ N]{U}e  dA  (4.15) 

A 

= {5  U)  ®T  [K]  e{U\e 


The  matrix  [K]e  is  the  element  elastic  stiffness  matrix.  A 
3*3  Gaussian  Quadrature  rule  can  be  used  to  reliably  evaluate 
the  element  stiffness  matrix.  The  3*3  integration  scheme 
evaluates  the  stiffness  matrix  exactly  for  parallelogram 
shaped  elements. 

When  thin  plates  are  modeled  using  this  element,  shear 
locking,  which  is  inherently  present  in  such  formulations  can 
cause  the  solution  to  be  very  stiff.  Reduced-order 
integration  is  often  used  to  avoid  this  problem,  but  the 
inexact  evaluation  of  the  stiffness  matrix  may  introduce 
spurious  modes  which  can  corrupt  the  solution.  The  problem  is 
well  studied  and  a fair  amount  of  work  is  currently  in 
progress  in  this  area  [57-60].  In  this  study,  no  further  work 
was  done  in  the  area;  however,  example  problems  which  show 
that  shear  locking  is  absent  for  the  geometries  that  are 
considered  in  this  study  are  discussed  in  Chapter  7.  The 
plate  element  discussed  here  is  implemented  in  UFPAC  such  that 
the  use  of  uniform  reduced  integration  is  optional.  Similar 
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isotropic  plate  elements,  however,  are  known  to  be  divergent 
in  the  thin  plate  limit  even  when  reduced  integration  is  used. 
Therefore,  care  should  be  exercized  when  using  this  element 
for  the  analysis  of  thin  plates. 

4.1.2  Initial  Stress  Stiffness  Matrix 

The  presence  of  an  initial  stress  in  a structure  modifies 
the  stiffness  of  the  structure.  Physically,  this  effect 
represents  the  coupling  between  the  inplane  and  out  of  plane 
deformations  and  is  modeled  as  a higher  order  effect. 

The  initial  stress  stiffness  matrix,  also  called  the 
geometric  stiffness  matrix  is  calculated  from  the  second  term 
in  the  incremental  virtual  work  expression.  The  second  term 
in  Equation  3.13  can  be  written  in  matrix  form  as 


(4.16) 


v. 
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where  the  vector  E is 
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(4.17) 


and  the  matrix  T is 


= {El  } -z  {E2  } 


[T]  = 


TOO 
0 T 0 
0 0 T 


(4.18) 


where  T is  the  initial  stress  matrix.  For  the  laminated  plate 
this  relation  can  be  expressed  in  terms  of  the  plate  force  and 
moment  resultants.  By  using  the  relations  in  Equations  4.17 
and  4.18,  Equation  4.16  takes  the  form 
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J{6£  }t  [T  ] {E}dv  = 
v0 

f{ 6F±  }t  [N  ] {£,}  dA  - f{8E2  }T  [M  ] dA  (4.19) 

A A 

- f { 8E [ } [M  ] {E2  } dA  + f {5E2  } [Z  ] { E2 } dA 

A A 

where 

Nx  0 0 

0 Nx  0 
0 0 Nx 


Mx  0 0 
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0 0 Mx 
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*••1 


[M]  = 


■ / 


z[T  ] dz  = 


z°'l 


[iV]  = 


/ [T] 


dz  = 


N-l,  Mx  and  and  Zx  are  the  matrices  of  stress  resultants,  and 
their  definition  is  identical  to  the  definition  of  [N] , [M] 
and  [Z],  The  only  difference  is  that  the  matrix  [T]  is 


replaced  by  [T] . 
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Equation  4.19  can  be  written  concisely  as 


f{8E  }t[T  ]{E  } dV 

Vo 


6 E1 ' 
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f < 
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M Z 
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-6E2 
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(4.20) 


The  matrix  of  stress  resultants  is  calculated  from  a static 
analysis  due  to  the  initial  load.  By  substituting  the 
expression  for  assumed  displacements  given  in  Equation  4.10, 
the  vector  {Ex  | -E2)T  can  be  written  as 


[P  ] IN] 


(4.21) 


where  [P]  is  a matrix  of  partial  differential  operators  and 
[N]  is  the  matrix  of  interpolation  functions  given  in  Equation 


4.11. 
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The  matrix  [P]  is 


-^-  0 0 0 0 

dx 
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By  substituting  Equation  4.21  in  Equation  4.20  the  geometric 
stiffness  matrix  can  be  expressed  as 


= /{8£}T[r  ]{£}  dv 

Vo  v0 


= f{bU\eT[N]T[P  }T 

A 


'N  M 
M Z 


[ P ] [N]{U)e  dA 


= {6  U)*t[Kg]{U} 


(4.22) 


The  geometric  stiffness  matrix  [KG]  can  be  evaluated  by  using 
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the  Gaussian  quadrature  scheme,  as  in  the  case  of  the  elastic 
stiffness  matrix. 

4.1.3  Mass  Matrix 

The  mass  matrix  is  calculated  in  a similar  manner  from 
the  third  term  in  the  incremental  virtual  work  expression. 
The  difference  between  the  present  derivation  and  the 
formulation  for  conventional  elements  is  due  to  the  offset 
factor.  The  acceleration  field  is  approximated  by  using  the 
same  interpolation  functions  that  were  used  to  approximate  the 
displacement  field.  Within  each  element,  the  acceleration 
field  is  assumed  to  be  of  the  form 

{u}  = [N]  {U}e  (4.23) 

where  [N]  is  the  matrix  of  interpolation  functions,  0i;  and 
{U}e  is  the  vector  of  elemental  acceleration  degrees  of 
freedom.  The  nonzero  entries  of  [N]  are  given  in  Equation 
4.11,  and  (pi  are  defined  in  Equation  4.12.  By  using  Equation 
4.23,  the  third  term  in  the  incremental  virtual  work  equality 
can  be  written  in  matrix  form  as 


/ 90^uiu1  dV  = Jp0{6 U\T  [U  } dV 
Vo  Vo 


(4.24) 


47 


The  matrices  (<5U)  and  {U}  in  Equation  4.24  are  defined  as 


= [A] 


[6 U ] 


where 
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Equation  4.24  can  now  be  written  as 


Jp0{6t7}r{/7}  dV 
V0 

= J p0{6U}r[^]  T [ A]t  [A]  [i7]{u}  dV 


(4.25) 


={6  U]  [ M]{U  } 
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The  consistent  element  mass  matrix  can  be  calculated  from  the 
expression 


[M]  =hfp0[N]T[A]  [AT]  dA  (4.26) 

A 


The  matrix  [A]  will  be  a diagonal  matrix  if  the  reference 
plane  is  not  offset.  The  entries  in  [A]  are  given  by 


1 . e . , 
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12 


where  z0  is  the  offset  and  h the  plate  thickness. 
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4.2  Variable  20-Node  Prestressed  Isoparametric  Solid 

The  three-dimensional , variable-twenty-node  isoparametric 
solid  is  an  isotropic  elastic  element  containing  the  linear- 
prestressed  nonlinearity  described  in  Chapter  3.  The  element 
is  suitable  for  both  static  and  dynamic  analyses  that  include 
this  type  of  nonlinearity.  Like  the  variable-node  plate 
element,  this  element  is  based  on  variable  order  interpolation 
functions  for  both  the  coordinates  defining  the  element  and 
the  unknown  (displacement)  variable.  The  element  may  contain 
a minimum  of  eight  nodes  to  a maximum  of  twenty  nodes,  with 
the  eight  node  element  being  the  basic  trilinear  solid.  Any 
of  the  nodes  from  9 to  20  can  be  added  to  introduce  a higher 
order  effect  (quadratic)  along  an  edge.  The  element  is 
called  a "serendipity  element"  since  it  has  no  interior  nodes, 
or,  nodes  that  are  not  shared  between  neighboring  elements. 

The  location  of  the  nodes  in  the  local  coordinate  system 
are  shown  in  Figure  4.2.  Since  each  of  the  nodes  from  9 
through  20  are  optional,  thousands  of  node  patterns  can  be 
generated  by  a single  element.  The  flexibility  is  extremely 
useful  in  making  transitions  from  coarse  to  fine  meshes,  and 
in  keeping  the  number  of  active  degrees  of  freedom  down  to  the 
minimum  required  even  when  uniform  meshes  are  used.  The 
element,  which  can  have  4 to  8 nodes  per  face,  is  fully 
compatible  with  the  plate  element  described  in  Section  4.1.1 
if  the  offset  parameter  for  the  plate  element  is  correctly 
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Figure  4.2.  Variable  Node  Solid  Element 
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defined.  Although,  the  three-dimensional  continuum  is  usually 
used  in  the  analysis  of  solids  and  thick  walled  structures,  it 
is  used  effectively  as  a shear  film  in  the  present  study  to 
model  the  shear  deformation  of  the  viscoelatic  damping  layer 
shown  in  Figure  2.1. 

4.2.1  Elastic  Stiffness  Matrix  and  Mass  Matrix 

The  variable  order  of  the  node  numbering  in  the  local 
coordinate  system  is  as  shown  in  Figure  4.2.  A concise 
definition  of  the  interpolation  functions  for  the  above  node 
numbering  is  given  in  Bathe  [44].  The  solid  element  is  used 
with  the  plate  element,  which  in  general,  has  six  degrees  of 
freedom  per  node  in  the  global  coordinate  system;  therefore, 
the  stiffness  matrix  of  the  20-node  solid  has  to  be  formulated 
accordingly.  A detailed  development  of  the  elastic  element 
stiffness  matrix  and  the  mass  matrix  can  be  found  in  Bathe 
[44],  and  only  a brief  derivation  is  included  in  this 
discussion. 

Within  each  element,  the  displacement  field  is  assumed  to 
be  of  the  form 

{u}  = [AT]  (4.27) 

where  [N]  is  the  matrix  of  interpolation  functions,  01 , and 
{ U } e is  the  vector  of  elemental  degrees  of  freedom.  The 
nonzero  entries  of  [N]  are 
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Nki  = 4>ic  (4.28) 
(k=l , 3 , (j=  1,20  ( i= ( j -1) *6  + 1 ))) 

[P]  is  the  matrix  of  differential  operators  and  [G]  the 
constitutive  matrix  for  an  isotropic  material.  [P]  and  [G] 
are  defined  as 
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The  matrix  [G]  is  symmetric,  and  elements  not  shown  are 
zeroes.  The  element  elastic  stiffness  matrix,  after  making 
the  appropriate  substitutions  in  the  first  term  of  the  virtual 
work  expression,  is  derived  in  terms  of  [N] , [P]  and  [G]  as 
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[K\  e = f [N]  T [P]  T [G]  [P]  [N]  dV 


(4.29) 


The  mass  matrix,  after  making  appropriate  substitutions  in  the 
third  term  of  the  virtual  work  expression,  is  derived  as 


The  stiffness  and  mass  matrices  are  calculated  by  using  an 
appropriate  Gauss  guadrature  rule. 

4.2.2  Initial  Stress  Stiffness  Matrix 

The  initial  stress  stiffness  matrix  of  the  solid  element 
is  calculated  from  the  second  term  in  the  incremental  virtual 
work  expression.  The  second  term  in  Equation  3.13  can  be 
written  in  matrix  form  as 


[M]e  = J p0  [N]  T [N]  dV 


(4.30) 


= |{£}r[T  ]{£■}  dV 


(4.31) 
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where  the  vector  E is 


and  the  matrix  T is 
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where  T is  the  initial  stress  matrix. 

By  making  the  necessary  substitutions,  the  geometric 
stiffness  matrix  of  the  solid  element  can  be  defined  as 


/ = /{5£}T[r  ]{£•}  dv 

vB 

= f{6U)eT[N]T[P  ]t[T][P  ] lN]{U}e  dA  (4.32) 

A 

= (6  U)er[Ka]{U) 

[P]  is  the  matrix  of  differential  operators  such  that 


{E  } = [P  ] [N  ] {U}e 


(4.33) 
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4.3  Laminated  Offset  Beam 


The  offset  beam  element  with  nodes  offset  to  the  top  is 
shown  in  Figure  4.3.  In  the  case  of  beams  treated  with 
constrained  layer  damping  treatments,  the  base  structure  and 
the  constraining  layer  were  modeled  by  using  the  three-node, 
seven-degree-of-freedom,  laminated  offset  beam  element.  The 
element  is  shear-deformable,  which  is  significant  when 
modeling  fiber-reinforced  composites.  A key  feature  of  this 
element,  like  the  plate  element  discussed  in  Section  4.1,  is 
its  ability  to  account  for  coupling  between  stretching  and 
bending  deformations.  The  offset  ability  allows  for  the  beam 
nodes  to  be  offset  to  one  surface  of  the  beam,  coincident  with 
the  nodes  of  the  adjoining  element  on  the  surface  of  the  beam. 

The  element  stiffness  matrix  of  the  offset  beam  element 
is  formulated  as  follows.  The  assumed  displacement  field  is 


u0  and  i|j  are  defined  by  using  linear  interpolation  functions, 


u(x,  z)  = u0  - z i|r  (x) 
w(x,  z)  = w{x) 
t (x,  z)  = i|r  (x) 


(4.34) 


u(x)  =[  ^ x't  0 
2 


(4 . 35) 


2 


0 


2 
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Figure  4.3.  Laminated  Offset  Beam  Element 
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where  u:,  u3,  1(1!  and  i|i3  are  corresponding  nodal  displacements, 
w is  defined  by  using  quadratic  interpolation  functions, 


' x(x-l) 

T 

' 

2 

Wi 

w = 

( 1-x 2) 

w2e- 

x(x+l) 

w® 

2 

where  w®,  w|  and  w3  are  nodal  displacements.  By  using  higher 
order  interpolation  functions  for  w,  shear  locking,  which  is 
inherently  present  in  shear  deformable  beam  elements,  can  be 
avoided. 

Strains  are  expressed  in  terms  of  displacements  as 


The  stress-strain  relation  is  of  the  form 


(4.37) 


(4 . 38) 


By  using  the  above  relations,  the  first  term  of  the  virtual 
work  expression  can  be  written  in  matrix  form  as 


f Cijxsers?>  eisdv  = j {bE)[G]{E)  dA  (4.39) 

Vo  A 


where  the  vector  { E ) and  the  matrix  [G]  are 
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{E\  = Wo.x  ~ ^ 


ABO 


[G\  = B D 0 
OOF 


z0+h/2 


A,  B,  D,  F = J Q,  z Q,  z2  Q,  Qs  dz 


za-h/2 


The  element  elastic  stiffness  matrix  [K]e  is  calculated  by 
substituting  the  expressions  for  the  assumed  displacement 
field  in  the  first  term  of  the  virtual  work  expression.  The 
element  stiffness  matrix  is  defined  such  that  the  first  term 
in  the  incremental  virtual  work  expression  reduces  to 


The  mass  matrix  and  geometric  stiffness  matrix  of  the 
laminated  beam  element  are  calculated  as  described  in  Section 
4.1  for  the  plate  element.  As  can  be  seen  from  the  derivation 
of  the  element  stiffness  matrix,  the  procedure  is  just  a 
specialization  of  the  derivation  for  the  plate  elements. 


(4.40) 


Vr 


A 


= {6ue}[Fe]{ue} 
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4.4  Plane  Solid  Element 

The  plane  solid  element  was  derived  to  model  the 
viscoelastic  damping  layer  in  damped  sandwich  beams; 
therefore,  it  is  required  to  be  compatible  with  the  laminated 
beam  element  of  Section  4.3.  The  plane  element  is  a six-node, 
ten-degree-of-freedom  element  and  is  shown  in  Figure  4.4.  In 
order  to  be  compatible  with  the  beam  element,  the  order  of 
interpolation  functions  was  chosen  as  described  below,  x- 
displacement  and  the  coordinates  were  interpolated  by  using 
linear  interpolation  functions.  The  nodal  displacements  in 
the  z direction,  w,  were  interpolated  using  quadratic 
interpolation  functions  in  the  x-direction  and  linear 
interpolation  functions  in  the  z-direction.  The  use  of  higher 
order  interpolation  for  w improves  the  performance  of  the 
element  in  situations  where  the  loading  causes  w to  be  much 
larger  than  u,  as  in  beam  bending  problems. 

As  for  the  other  elements  discussed  in  this  chapter, 
the  element  matrices  are  derived  by  appropriate  discretization 
of  the  virtual  work  expression.  The  derivation  of  the  element 
matrices  is  simply  a two-dimensional  version  of  the  procedure 
presented  in  Section  4.2  for  the  three-dimensional  solid 
element . 

By  using  the  elements  formulated  in  this  chapter,  beams 
and  plates  treated  with  a constrained  viscoelastic  layer  can 
be  modeled  as  shown  in  Figures  4.5  and  4.6.  In  each  case,  the 
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Figure  4.4 


Two-dimensional  Plane  Solid 
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Figure  4.5. 


Damped  Beam  Finite  Element  Model 
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Figure  4.6.  Damped  Plate  Finite  Element  Model 
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three  layer  sandwich  can  be  modeled  by  just  two  layers  of 
nodes.  In  general,  as  shown  in  Figure  4.6,  both  the  outer 
layers  can  be  arbitrarily  laminated. 


CHAPTER  5 

CALCULATION  OF  LOSS  FACTOR 

Damping  is  an  important  consideration  in  several 
applications,  and  better  results  are  often  achieved  from 
damping  by  design.  For  low  levels  of  damping,  the  choice  of 
method  is  influenced  more  by  efficiency  than  by  difference  in 
numerical  result.  In  this  chapter,  the  mathematical  model 
used  to  account  for  damping  and  two  different  methods  to 
evaluate  system  loss  factor  are  discussed.  In  both  methods 
the  overall  damping  of  the  structure  is  calculated  from  the 
independent  properties  of  the  constituent  materials  that  form 
the  composite  structure.  An  efficient  modification  of  the 
modal  strain  energy  method  for  the  estimation  of  damping 
ability  about  a prestressed  configuration  is  also  presented. 

5.1  Complex  Modulus  Approach 
There  are  many  guantitative  measures  of  structural 
damping  ability  based  on  the  effect  it  has  on  the  response  of 
the  system.  Most  of  the  measures  of  damping  that  are  used 
currently,  such  as  resonance  bluntness,  half-power  bandwidth, 
logarithmic  decrement,  percentage  of  critical  damping  etc,  are 
defined  based  on  linear  single  degree  of  freedom  systems  with 
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frequency  independent  parameters.  Of  the  many  different 
measures  of  damping  used  in  resonant  analysis,  loss  factor  and 
quantities  closely  related  to  it  are  the  most  generally 
applicable.  Loss  factor  may  be  defined  in  terms  of  energy 
quantities  so  as  to  apply  to  nonlinear  systems  and  systems 
with  frequency-dependent  parameters.  All  the  different 
measures  of  damping,  under  the  assumption  of  low  damping,  are 
used  interchangeably  and  are  related  as  shown  in  Table  5.1 
[61]  . 

The  concept  of  structural  damping,  also  called  hysteretic 
damping,  was  originally  postulated  as  a basis  for  describing 
the  internal  damping  properties  of  solids  [62,63].  Mycklestad 
[64]  and  others  [65-67]  showed  that  if  the  hysteresis  loop  for 
a material  has  a certain  shape,  then  the  damping  can  be 
considered  adequately  by  using  a complex  stiffness,  K’  + iK  , 
where  the  real  part  is  called  the  storage  modulus  and  the 
imaginary  part  the  loss  modulus.  Loss  factor,  a measure  of 
damping,  is  the  ratio  of  the  imaginary  part  to  the  real  part. 
The  complex  modulus  is  only  defined  for  harmonic  excitation; 
however,  since  any  real  signal  can  be  represented  as  a Fourier 
sum,  the  complex  modulus  approach  can  also  be  used  for  the 
analysis  of  a system  subjected  to  arbitrary  excitation. 

From  a rigorous  mathematical  view  point,  hysteretic 
damping  is  not  a physically  realizable  energy  dissipation 
mechanism.  It  is  known  that  complex  stiffnesses  which  do  not 
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Table  5.1.  Relation  Between  Damping  Measures 


Measure 

Damping 

Ratio 

Loss 

Factor 

Log  De- 
crement 

Quality 

Factor 

Speci- 

fic 

Damping 

Damping 

Ratio 

C 

77/2 

A/7 t 

1/2Q 

D/47TW 

Loss 

Factor 

2 C 

rj 

2A/7T 

l/Q 

D/27TW 

Log  De- 
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7rC 

ITT] / 2 

A 

7T/2Q 

D/4W 

Quality 
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1/2  C 

1/r? 

7T/2  A 

Q 

27TW/D 
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D 
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vary  with  frequency  lead  to  noncausal  system  response 
(response  anticipates  input)  [68],  and  a finite  energy 
dissipation  rate  at  zero  frequency.  This  matter  has  been 
stronqly  debated  for  several  years;  however,  if  the  measured 
physical  properties  are  used,  the  problem  disappears. 

The  loss  factor,  77,  of  a structure  executing  steady  state 
vibrations  is  defined  in  terms  of  energy  quantities  as 


D is  the  energy  dissipated  per  cycle  or  the  energy  that  has  to 
be  supplied  to  sustain  the  steady  state  vibration.  W is  the 
maximum  energy  associated  with  vibration.  Ungar  and  Kerwin 
[69]  correctly  observed  that  the  definition  of  W was 
unambiguous  only  for  a lightly  damped  system.  They  defined  W 
as  the  total  energy  stored  in  an  identical  lossless  system. 
They  also  showed  that  for  an  arbitrary  series-parallel  network 
of  viscoelastic  springs,  the  composite  loss  factor  reduced  to 
the  simple  expression 


n = 


E K 

EWn 


(5.2) 


where,  rjn  and  Wn  are  the  loss  factors  and  the  maximum  strain 
energy  stored  in  the  nth  spring  respectively. 


Hysteretic  damping  in  a continuous  composite  system  can 
be  modeled  by  using  the  correspondence  principle  of 
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viscoelasticity  [70].  According  to  this  principle,  the 
elastic  constants  of  the  materials  are  simply  replaced  by  the 
corresponding  viscoelastic  counterparts.  Two  approaches  to 
calculating  loss  factors  and  resonant  frequencies  from  the 
resulting  complex  valued  system  of  equations,  in  the  context 
of  implementation  by  the  finite  element  method,  are  discussed 
in  the  following  sections. 

5.2  Techniques  for  the  Estimation  of  Loss  Factor 
The  two  approaches  considered  are  called  the  direct 
frequency  response  method  and  the  modal  strain  energy  method 
respectively.  The  methods  and  their  relative  benefits  are 
discussed  in  the  remainder  of  this  chapter. 

5.2.1  Direct  Frequency  Response 

In  the  direct  frequency  response  method,  a forced 
vibration  over  a range  of  frequencies  is  considered,  and  a 
theoretical  response  spectrum  is  generated  from  which  the 
overall  damping  of  the  system  is  calculated.  For  the  forced 
vibration  problem,  the  objective  is  to  predict  the  linear, 
damped,  steady-state  response  of  structures  about  a linear 
preloaded  equilibrium  configuration.  Since  the  complex 
modulus  approach  is  used,  the  governing  differential  equations 
are  complex.  Following  finite  element  discretization,  the 
equations  governing  the  time  dependent  response  of  a 
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prestressed  body  are  obtained  as 

[ [ K ] + [Ka]  ] {iJl  + [M] 


(5.3) 


where 

[K]  = global  elastic  stiffness  matrix  (complex) 

[Kg]  = global  geometric  stiffness  matrix 

[M]  = global  mass  matrix 

{ U } = global  displacement  vector  (complex) 

(U)  = global  acceleration  vector  (complex) 

A harmonic  excitation  force  is  considered,  i.e., 

F = feiot  (5.4) 

where,  w is  the  forcing  frequency,  t is  the  time  and  i=(-l)1/2. 
The  response  due  to  the  applied  force  is  assumed  to  be 
harmonic  and  vibrating  at  the  excitation  frequency.  The 
equation  of  motion  reduces  to 

[ [K\  + [Kg]  ] (d-w2[W]  fo}  = {f}  (5.5) 

The  matrix  on  the  left  hand  side  is  called  the  displacement 
impedance  matrix.  For  each  frequency  of  excitation,  the 
displacement  per  unit  applied  force  is  calculated  by  solving 
the  system  of  complex-valued  simultaneous  linear  equations, 
and  the  response  function  is  thereby  generated.  Loss  factor 
is  calculated  from  the  generated  response  spectrum  by  using 
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the  half-power-bandwidth  technique,  or  from  the  real  part  of 
the  spectrum  as  shown  in  Figures  5.1  and  5.2  respectively. 

Loss  factor  can  also  be  calculated  by  using  the  energy 
ratio  method  described  in  the  previous  section.  After  the 
resonant  frequency  and  the  corresponding  displaced  shape  at 
resonance  have  been  calculated,  the  loss  factor  can  be 
calculated  from 


u'tk'u 

where  U and  u*  are  the  complex  displacement  and  the  complex 
conjugate  of  the  displacement  respectively. 

The  direct  frequency  response  has  many  drawbacks.  It  is 
extremely  expensive  because  the  displacement  impedance  matrix 
has  to  be  formed  at  each  frequency  for  which  displacements  are 
desired.  Refined  frequency  resolution  is  required  and  this 
poses  problems  such  as  recalculation  of  the  stiffness  matrix 
(the  viscoelastic  material  has  frequency  dependent 
properties) . If  the  system  is  highly  damped,  coupling  between 
modes  may  be  present,  and  therefore  it  may  not  be  possible  to 
excite  one  mode  without  exciting  the  other.  In  this  case  the 
deflected  shape  will  not  necessarily  resemble  the  mode  shape, 
and  the  calculated  damping  is  not  a true  modal  value. 
Finally,  the  information  provided  by  the  direct  frequency 
response  method  is  not  useful  to  a designer  in  improving  the 
performance  of  a candidate  design.  The  technique,  however, 
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Figure  5.1.  Real  Part  of  the  Response  Spectrum 
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Figure  5.2.  Half  Power  Bandwidth  Method 
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serves  to  verify  the  predictions  of  the  other  methods  since  it 
makes  no  assumptions  about  the  level  of  damping. 

5.2.2  Modal  Strain  Energy  Approach 

The  second  technique  is  based  on  an  eigenvalue  problem 
which  derives  from  setting  the  forcing  terms  to  zero,  i.e., 

[ IK]  + [Kg]  ] (d}  = G)2  [M]  {l/t  (5.7) 

The  analysis  can  be  performed  based  on  the  associated  complex 
eigenvalue  problem  or  the  simplified  real  eigenvalue  problem. 
In  the  complex  eigenvalue  method,  the  resonant  frequencies  and 
complex  mode  shapes  of  the  damped  system  are  determined  by 
solving  the  complex  eigenvalue  problem.  The  modal  loss  factor 
can  be  calculated  from  Equation  5.6,  but  in  this  case,  U is 
the  complex  mode  shape.  The  complex  eigenvalue  method  is 
numerically  inefficient  and  is  usually  not  used. 

In  the  method  based  on  the  real  eigenvalue  solution, 
natural  frequencies  and  mode  shapes  are  calculated  for  the 
undamped  system.  The  stiffness  matrix  and  the  corresponding 
nodal  displacements  are  real  since  only  the  real  eigenvalue 
problem  is  solved.  This  technique  is  valid  only  for  systems 
with  relatively  small  levels  of  damping  where  the  mode  shapes 
and  frequencies  of  the  damped  and  the  undamped  structure  are 
similar.  This  causes  the  fraction  of  strain  energy  stored  in 
each  element  to  also  be  similar.  The  system  loss  factor  is 
calculated  as  the  weighted  sum  of  the  loss  factors  of  the 
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individual  elements,  where  the  weighting  factor  is  the 
fraction  of  the  strain  energy  stored  in  the  element. 

If  the  prestress  present  in  a structure  causes  a 
stiffness  increase,  then  the  energy  stored  per  cycle  can 
change  by  a significant  amount.  The  weighting  factors  have  to 
be  modified  to  account  for  this  increase  in  stored  energy, 
since  the  dissipated  energy  per  cycle  does  not  change  very 
much.  The  loss  factor  for  vibration  about  a prestressed 
configuration  is  calculated  by 

r,  = niSf/SiS  Et  (5.8) 

where 

n = total  number  of  elements 

Vi  = loss  factor  of  the  ith  element 

E?  = elastic  strain  energy  stored  in  ith  element 

calculated  as 

1/2  UT  [Jr]  (u) 

Ei  = sum  of  elastic  strain  energy  and  strain  energy 
due  to  preload  in  the  ith  element  calculated  as 

1/2  (uK  {u} 

The  modal  strain  energy  based  method  is  the  most  popular 
because  of  the  computational  efficiency  of  the  method.  A 
variation  of  the  modal  strain  energy  method  has  been  used  by 
Soni  [71]  to  calculate  loss  factors  of  prestressed  structures. 
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In  his  approach,  the  resonant  frequency  is  first  located  by 
using  direct  frequency  response  as  described  in  Section  5.2.1. 
The  deflected  shape  of  the  damped  structure  at  resonance  is 
then  used  to  calculate  the  strain  energy  fractions  stored  in 
the  element.  If  this  approach  is  used,  the  computational 
benefits  of  the  modal  strain  energy  method  are  lost  due  to 
reasons  discussed  under  the  direct  frequency  response  method. 

Even  though  the  assumption  of  small  damping  is  made  in 
the  modal  strain  energy  approach,  accurate  results  have  been 
reported  for  core  loss  factors  as  high  as  1 [38].  However, 
for  higher  core  loss  factors  the  modal  values  calculated  by 
the  modal  strain  energy  method  will  depart  from  those 
calculated  by  the  other  methods,  and  the  technique  chosen  must 
be  evaluated  carefully.  Also,  in  the  case  of  complex 
geometries,  the  vibration  mode  shape  at  resonance  will  not 
necessarily  correspond  to  the  deflected  shape  due  to  a point 
forcing  function.  This  results  in  larger  differences  in  the 
predicted  loss  factors  between  the  two  methods.  In  general, 
the  modal  methods  provide  information  about  a basic  and  more 
general  behavior  of  the  damped  structure,  over  a wide 
frequency  range,  without  accounting  for  the  actual  loading 
situation.  The  direct  frequency  response  method  will  provide 
detailed  response  information  in  the  specified  range  of 


interest . 
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The  variation  of  material  properties  of  the  viscoelastic 
damping  material  with  frequency  is  very  pronounced,  and  in 
eigenvalue-based  techniques,  the  frequencies  of  interest  are 
not  known  a priori.  Therefore,  an  iterative  procedure  has  to 
be  used,  starting  with  a rough  estimate  of  resonant 
frequencies.  The  material  properties  of  the  viscoelastic 
material  are  updated  corresponding  to  the  calculated  resonant 
frequencies.  The  procedure  is  repeated  to  convergence  for 
each  mode  of  interest.  System  loss  factors  calculated  without 
iterating  will  produce  correspondingly  inaccurate  results. 
This  is  one  problem  which  is  avoided  in  the  direct  frequency 
method  since  the  loading  frequency  is  known  ahead  of  time. 


CHAPTER  6 

EXPERIMENTAL  MEASUREMENT  OF  LOSS  FACTOR 

The  most  common  methods  used  to  measure  damping  are  the 
free  vibration  decay  method,  the  resonant  dwell  method,  the 
hysteresis  loop  method  and  the  frequency-response  technique 
[72-80].  For  all  experimental  measurements  of  loss  factor 
presented  in  this  paper,  the  frequency-response  technique[76- 
80]  was  used.  This  technique  offers  potential  for  rapid  non- 
destructive evaluation  of  materials  and  structures. 

In  this  technique,  the  specimen  is  excited  impulsively 
with  a controlled-impact  hammer  which  has  a force  transducer 
attached  to  its  head.  The  specimen  response  is  sensed  by  a 
non-contacting  eddy  current  proximity  probe.  The  signals  from 
the  force  transducer  and  the  motion  transducer  are  input  to  a 
Fast  Fourier  Transform  ( FFT)  analyzer  which  displays  the 
frequency  response  spectrum.  A block  diagram  of  the 
instrumentation  required  for  the  flexural  vibration  test  is 
shown  in  Figure  6.1.  By  analyzing  the  resonant  peak  for  a 
particular  mode,  the  loss  factor,  a measure  of  damping,  is 
obtained  from  the  real  part  of  the  response  spectrum  as 
explained  in  Figure  5.1.  Loss  factor  may  also  be  calculated 
by  using  the  half-power  bandwidth  technique  which  is  shown  in 
Figure  5.2. 
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Figure  6.1.  Flexural  Vibration  Test  Apparatus 
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One  feature  of  this  technique  is  that  the  excitation 
level  is  accurately  controlled,  therefore,  the  amplitude  of 
vibration  of  the  specimen  can  be  kept  to  a minimum  (thereby 
reducing  air  damping  to  a minimum) . Also,  the  response 
function,  which  is  identical  in  shape  to  the  transfer  function 
after  ensemble  averaging  can  be  used  for  damping  measurements 
[80]  . This  result  is  valid  for  impulsive  excitation  only. 
Details  about  the  experimental  method  and  other  techniques  for 
measuring  damping  are  available  in  the  papers  referenced  in 
this  section.  The  technique  was  used  mainly  to  generate 
material  property  data  which  was  required  as  input  to  the 
finite  element  model.  Some  comparisons  between  experimental 
and  analytical  predictions  for  loss  factors  of  beams  with 
constrained  viscoelastic  damping  treatment  applied  to  its 
surface  were  also  made. 


CHAPTER  7 

VERIFICATION  PROBLEMS 

7 . 1 Static  Displacement  of  a Laminated  Beam 
The  stiffness  matrix  formulation  of  the  beam  element  is 
verifi-ed  in  this  section  by  examining  the  performance  of  the 
element  in  different  situations.  The  finite  element 
calculations  of  beam  deflection  are  compared  with  the  Euler- 
Bernoulli  solution  for  different  beam  thicknesses  to  show  that 
the  element  performs  well  in  the  thin  beam  limit.  That  is,  as 
the  thickness  of  the  beam  approaches  zero,  the  finite  element 
predictions  approach  the  predictions  of  the  Euler-Bernoulli 
beam  theory. 

Figure  7.1  shows  that  there  is  good  agreement  between  the 
Euler-Bernoulli  theory  solution  for  large  aspect  ratios  (beam 
length/thickness) . As  the  aspect  ratio  decreases,  the 
deflections  predicted  by  the  shear-deformable  finite  element 
are  considerably  larger  than  the  Euler-Bernoulli  predictions. 
Generally,  it  is  very  important  to  consider  shear  deformation 
in  composite  beams  because  of  their  low  shear  moduli.  The 
finite  element  solutions  were  also  compared  with  Euler- 
Bernoulli  solutions  for  beams  with  "almost  zero"  thicknesses 
and  were  found  to  be  exact.  Shear  locking,  which  is 
inherently  present  in  such  formulations  is  completely  avoided 
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Figure  7.1.  Static  Displacement  Vs  Aspect  Ratio  (Beam) 
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because  a higher  order  interpolation  function  is  used  for 
deflection  ( z-displacement) . 

The  working  of  the  offset  capability  of  the  element  was 
verified  by  comparing  the  solutions  from  a model  with  one 
element  through  thickness,  and,  several  models  constructed  as 
a sandwich  of  top  and  bottom  offset  beams.  The  results  were 
found  to  be  identical  as  should  be  the  case.  The  beam  in  the 
current  problem  is  fixed  at  one  end,  and  its  properties  are 
E = 127.9  GPa 
G = 7.6  GPa 

Length  = 0.01  to  0.2  m (ten  equal  length  elements) 

Width  = 0.01  m 
Thickness=  0.01  m 

7.2  Static  Displacement  of  a Laminated  Plate 

In  this  section,  finite  element  solutions  are  compared 
with  analytical  solutions  for  several  simply  supported 
laminated  plates  subjected  to  a concentrated  center  load.  The 
analytical  solution  is  a fifty-term  series  solution  based  on 
a shear-deformable  plate  theory.  The  theory  is  an  extension 
of  the  one  developed  by  Reissner  [81],  and  Mindlin  [82]  for 
homogenous  isotropic  plates.  The  extension  was  originally  due 
to  Yang  et  al.  [83],  and  was  modified  later  by  Whitney  and 
Pagano  [84].  The  governing  differential  equations  for  bending 
of  symmetric  laminates  under  lateral  loads  can  be  found  in 
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Whitney  and  Pagano  [84].  A comparison  between  the  results 
from  the  present  analysis  and  the  series  solution  is  given  in 
Table  7.1. 

The  plate  properties  and  dimensions  are 
Ei  = 14  0 GPa  E2  = 15  GPa  G12  = 6 GPa 
v 12  = 0.21  v23  = 0.21 

Ply  thickness  = 0.0005  m 

Four  equal  sized  elements  were  used  along  each  edge  to  model 
the  one  quarter  of  the  plate.  For  the  first  set  of  results 
the  aspect  ratio  (length  or  width  divided  by  the  thickness)  is 
8.33  and  in  the  second  case  the  aspect  ratio  is  83.33.  Finite 
element  calculations  for  each  set  of  results  were  made  for  two 
cases.  In  the  first  case,  a 3X3  integration  scheme  was  used 
to  evaluate  the  stiffness  matrix,  and  in  the  second,  a uniform 
reduced  scheme  (2X2)  was  used  to  calculate  the  element  matrix. 
Percent  error  as  compared  to  the  series  solution,  when  a (3X3) 
scheme  is  used  is  given  in  Table  7.1.  Finite  element 
solutions  are  almost  identical  to  the  series  solution  when 
uniform  reduced  integration  is  used  to  calculate  the  stiffness 


matrix. 
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Table  7.1.  Static  Displacements  of  the  Laminated  Plates 


(length  * width  * thickness)  = 

0.2  * 0.2  * 0.003;  Aspect  Ratio  = 8.3 

Stacking 

Sequence 

Series 

Solution 

Integration 
Order  (3x3) 

Integration 
Order  (2  x 2) 

[0]6 

1.908  E-5 

1.895  E-5 
(3.8%) 

1.904  E-5 

[0/902]s 

1.756  E-5 

1.719  E-5 
(2 . 1%) 

1.755  E-5 

ID 

0 
cr> 

1  i 

1.908  E-5 

1.835  E-5 
(3.8%) 

1.904  E-5 

[ 902/  0 ] s 

1.862  E-5 

1.801  E-5 
(3 . 3%) 

1.860  E-5 

(length  * width  * thickness)  = 

2 * 2 * 0.003;  Aspect  Ratio  = 83.33 

Stacking 

Sequence 

Series 

Solution 

3x3 

2x2 

[0]6 

1.880  E-3 

1.711  E-3 

(9%) 

1.851  E-3 

[ 0/9  02  ] s 

1.733  E-3 

1.626  E-3 
(6.2%) 

1.719  E-3 

C 90  ] 6 

1.880  E-3 

1.711  E-3 
(9%) 

1.851  E-3 

C 9 02/  0 ] s 

1.835  E-3 

1.686  E-3 
(8%) 

1.811  E-3 
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7.3  Loss  Factors  and  Frequencies  of  a Sandwich  Beam 
7.3.1  Sandwich  Beam  with  No  Prestress 

The  resonant  frequencies  and  loss  factors  of  a sandwich 
beam,  calculated  by  using  the  finite  element  model  developed 
in  this  study,  were  compared  with  the  results  obtained  in 
previous  studies.  This  problem  has  been  studied  extensively 
in  other  works  [39,40].  The  calculated  resonant  frequencies 
and  beam  loss  factors  are  compared  with  the  finite  element 
solution  of  Soni  and  Bogner  [39]  and  a sixth  order  beam  theory 
solution  [26].  The  geometry  of  the  beam  is  shown  in  Figure 
7.2.  For  the  purpose  of  comparison,  a hypothetical  damping 
material  with  frequency  and  temperature  invariant  properties, 
similar  to  the  one  considered  by  Soni  and  Bogner  [39],  is 
used.  The  material  properties  of  the  sandwich  beam  are  given 
in  Table  7.2. 

The  damped  frequencies  and  corresponding  modal  loss 
factors  from  the  present  finite  element  analysis  are  compared 
with  those  from  previous  analyses  in  Table  7.3.  The  values  of 
the  damped  resonant  frequencies  and  the  system  loss  factors 
are  in  good  agreement.  Loss  factor  was  calculated  by  using 
the  Direct  Frequency  Response  Method  described  in  Chapter  5. 
In  the  present  model,  a total  of  80  active  degrees  of  freedom 
was  used  to  model  the  sandwich  beam.  Soni  and  Bogner' s model 
employed  440  active  degrees  of  freedom. 


Figure  7.2.  Verification  Problem  7.3.1 
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Table  7.2.  Material  Properties  for  the  Sandwich  Beam 


ALUMINUM 

VISCOELASTIC 

CORE 

Thickness  (mm) 

1.524 

0.127 

Young's  modulus 
(GPa) 

69 . 000 

0.0021 

Poisson  ratio 

0.300 

0.499 

Loss  factor 

0.000 

1.000 

Specific 

gravity 

2 .800 

0.970 

88 


Table  7.3.  Results  (Section  7.3.1) 


Mode 

Loss  factor 

Frequency  (Hz) 

Present 

6th 

order 

Soni 

Present 

6th 

order 

Soni 

1 

0.0281 

0 . 0282 

0.0282 

63.86 

64 . 08 

64 .20 

2 

0 . 0244 

0.0242 

0 . 0243 

296.1 

296.4 

297 . 0 

3 

0.0154 

0.0154 

0.0153 

745.2 

743.7 

747 . 2 

4 

0.0088 

0.0089 

0.0088 

1404 

1393 

1408 

5 

0.0056 

0.0057 

0.0056 

2295 

2261 

2304 

6 

0 . 0037 

0.0039 

0 . 0038 

3427 

3343 

3446 
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7.3.2  Sandwich  Beam  With  Prestress 

The  prestress  formulation  and  the  technique  for 
calculation  of  modal  loss  factor  about  prestressed 
configurations  are  verified  in  this  section.  Results  from  the 
present  analysis  for  loss  factor  as  a function  of  preload,  for 
the  sandwich  beam  of  Section  7.3.1,  are  compared  with  results 
form  a three-dimensional  model.  The  beam  is  clamped  at  one 
end,  and  has  an  axial  preload  acting  at  the  other  as  shown  in 
Figure  7.3.  Loss  factor  and  resonant  frequency  were 
calculated  by  using  the  two  methods  discussed  in  Chapter  5 and 
using  the  elements  derived  in  the  current  study,  and  from  a 
fully  three-dimensional  model  using  ANSYS[85].  For  the 
purpose  of  comparison  the  viscoelastic  material  with 
frequency-independent  properties  considered  in  the  previous 
example  was  used  in  this  problem.  Resonant  frequencies  and 
loss  factors  are  compared  for  the  first  two  modes. 

The  properties  and  dimensions  of  the  beam  are  given  in 
Table  7.2  and  Figure  7.3.  A low  damping  core  was  used  since 
modal  strain  energy  method  is  accurate  when  core  loss  factors 
are  low.  The  results  are  listed  in  Table  7.3  and  show  that 
for  simple  geometries,  identical  loss  factors  should  be 
expected  from  the  two  methods,  provided  the  core  damping  is 
low.  For  core  loss  factors  greater  than  one  the  predictions 
of  the  modal  strain  energy  method  were  found  to  be 
consistently  higher  for  the  first  few  modes  of  vibration. 


\\\\\\ 


90 


Length  = 0.1778  m 


F Sin  (wt) 


Width  = 0.0254  m 
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Figure  7.3. 


Verification  Problem  7.3.2 
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Table  7.4.  Results  (Section  7.3.2) 


ANSYS  (3D) 

M.S.E. 

D.F.R. 

Mode 

PRELOAD 

(N) 

FREQ 

(Hz) 

LOSS 

FACTOR 

FREQ 

(Hz) 

LOSS 

FACTOR 

LOSS 

FACTOR 

1 

10 

66 

0.026 

66 

0 . 026 

0.026 

2 

303 

0.023 

303 

0.023 

0.024 

1 

100 

83 

0.017 

83 

0.017 

0.017 

2 

336 

0.018 

338 

0.018 

0.019 

1 

200 

98 

0 . 013 

97 

0.013 

0.013 

2 

375 

0.015 

375 

0.015 

0.015 

1 

300 

109 

0.011 

108 

0 . Oil 

0.011 

2 

403 

0.013 

401 

0.013 

0.013 
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7.4  Frequencies  and  Mode  Shapes  of  a Laminated  Plate 

In  this  section,  the  finite  element  predictions  of 
resonant  frequencies  and  mode  shapes  of  a 45  degree  plate  are 
compared  with  the  experimental  results  reported  by  Ashton  and 
Anderson  [86].  The  solutions  are  based  on  the  following  plate 
properties . 

E:  = 31  Msi  E2  = 2.7  Msi  G12  = 0.75  Msi 

v12  = 0.21  v23  = 0.21 

length  = width  = 12  in  thickness  = 0.424  in 
p = 1.92e-4  lb  secVin* 

The  results  for  the  first  five  vibration  modes  of  the 
plate  are  compared  in  Figures  7.4  through  7.8.  The 
experimentally  obtained  Chldani  figures  illustrate  node  lines 
of  the  characteristic  shapes.  Since  the  plate  is  strongly 
anisotropic,  the  node  lines  are  skewed  and  curved.  Careful 
examination  of  the  finite  element  results  show  that  the  mode 
shapes  are  reproduced  accurately  in  the  finite  element 
solution.  Ten  equal-size  elements  along  each  edge  were  used 
to  model  the  complete  plate.  The  agreement  between 
frequencies  is  reasonable,  and  the  difference  between  the 
frequencies  is  believed  to  be  due  to  the  edges  not  being 
perfectly  clamped  in  the  experimental  measurements. 

The  prestress  formulation  and  the  loss  factor 
calculations  for  the  three-dimensional  model  were  also 
verified  by  solving  the  sandwich  beam  problem  described  in 
Section  7.3.  Plate  and  three-dimensional  solid  finite 
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110.3  Hz 


Figure  7.4.  Laminated  Plate  (Mode  1) 
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186  Hz 


156  Hz 


Figure  7.5.  Laminated  Plate  (Mode  2) 
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256  Hz 


215  Hz 


Figure  7.6.  Laminated  Plate  (Mode  3) 
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280  Hz 


239  Hz 


Figure  7.7.  Laminated  Plate  (Mode  4) 
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380  HZ 


329  Hz 


Figure  7.8.  Laminated  Plate  (Mode  5) 
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elements  were  used  to  model  the  beam  as  shown  in  Figure  4.6. 
The  solutions  corresponding  to  the  bending  modes  were  found  to 
be  identical  to  the  predictions  of  the  two-dimensional  model 
given  in  Tables  7.2  and  7.3. 


CHAPTER  8 
EXAMPLE  PROBLEMS 

8.1  Loss  Factors  and  Frequencies  of  Some  Laminated  Composite 
Beams  with  Constrained-Layer  Damping  Treatments 

Two  different  materials  are  considered  for  the  purpose  of 
analysis.  Material  1 is  a hypothetical  material  with  low 
modulus  and  high  loss  factor,  typical  of  off-axis  composites. 
Typical  Glass-Epoxy  data  was  used  for  material  2.  The 
material  properties  of  material  1,  Glass-Epoxy  and  the  soft 
aluminum  constraining  layer  are  given  in  Table  8.1.  The 
damping  material  used  was  3M's  SJ2052x,  a class  of  constrained 
viscoelastic  damping  tape.  The  shear  modulus  and  loss  factor 
of  the  damping  material,  as  a function  of  temperature  and 
frequency  is  given  in  Figure  2.3. 

Structural  damping  with  and  without  the  add-on 
viscoelastic  layer  damping  was  evaluated  analytically  and 
experimentally  for  material  2.  Results  of  the  effects  of 
different  parameters  such  as,  the  quantity  of  treatment,  and 
location  of  treatment,  on  the  overall  damping  of  the  system, 
are  presented.  Three  different  unidirectional  composite 
specimens  were  tested,  and  the  average  measured  loss  factor 
was  used  as  input  data  to  the  finite  element  model.  The  same 
specimens  were  also  tested  after  application  of  the  damping 
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Table  8.1.  Material  Properties 


Material  1 Material  2 Constraining 

layer 

V,(%) 

50 

50 

P (g/ cm3) 

1.90 

1.90 

2.76 

El  (GPa) 

28 

36 

69 

Et  (GPa) 

8.8 

8 .80 

69 

^lt  (GPa) 

3 

3 

26 

VLT 

0.28 

0.28 

0.32 

0.01 

0.004 

0.005 

Hi 

0.015 

0.01 

0.005 
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treatment.  Each  specimen  was  203  mm  long  and  3.6  mm  thick. 
The  beam  configuration  is  shown  in  Figure  8.1. 

Figure  8.2  shows  the  magnitude  of  the  response  as  a 
function  of  the  frequency  of  the  forced  vibration  for  three 
different  lengths  of  treatment.  The  change  in  system  response 
due  to  addition  of  the  viscoelastic  material  can  be  seen  from 
the  figure.  Displacement  is  plotted  in  meters,  per  Newton  of 
applied  force.  Large  reductions  in  response  amplitude  can  be 
seen  due  to  application  of  the  damping  tape. 

Figure  8.3  shows  the  variation  of  loss  factor  with  tape 
length  for  material  1.  For  mode  1,  the  loss  factor  increases 
rapidly  from  b/L=0  to  b/L=0.6,  after  which  it  shows  a slight 
drop.  The  existence  of  a tape  length,  b,  for  which  b/L<l  and 
damping  is  optimal  is  significant.  This  result  shows  that 
shear  deformation  of  the  viscoelastic  core  is  the  primary 
source  for  energy  dissipation.  For  lengths  greater  than  the 
optimal  value,  the  deformation  of  the  viscoelastic  core 
follows  the  extensional  deformation  of  the  surface  of  the 
beam.  The  trend  observed  for  mode  2 is  different  from  that 
for  mode  1.  While  treatment  closer  to  the  root  of  the  beam 
seems  to  have  the  greatest  effect  on  mode  1,  the  center  of  the 
beam  seems  to  be  the  optimal  location  for  mode  2 . This 
suggests  that  application  of  damping  material  to  the  points  of 
high  strain  energy  density  yields  best  results.  These  results 
are  also  evident  from  Figures  8.4  and  8.5  which  show  the 
variation  of  system  loss  factor  with  tape  length,  for 
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Figure  8.1.  Beam  Configuration 


TIP  DISPLACEMENT  (m/N) 
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Figure  8.2.  Response  Spectra  for  Different  Lengths  of  Treatment. 
Material  1,  a/1  = 0 


LOSS  FACTOR 
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TAPE  LENGTH/BEAM  LENGTH 


Figure  8.3. 


Variation  of  Loss  Factor.  Material  1,  a/L  = 0 


LOSS  FACTOR 


TAPE  LENGTH/BEAM  LENGTH 


Figure  8.4.  Effects  of  Location  on  Loss  Factor 
(material  1,  mode  1) 


LOSS  FACTOR 
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TAPE  LENGTH/BEAM  LENGTH 


Figure  8.5.  Effects  of  Location  on  Loss  Factor 
(material  1,  mode  2) 
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different  locations  of  the  damping  treatment  for  mode  1 and 
mode  2 respectively. 

Figure  8.6  shows  the  analytical  and  experimental  results 
for  the  [ 0 ] 16  Glass-Epoxy  laminate  (Material  2).  Good 
agreement  between  the  data  is  seen  for  both  mode  1 and  mode  2. 
Analytical  and  experimental  results  project  similar  behavior 
for  varying  amounts  of  damping  treatment.  The  loss  factor  of 
the  base  structure  is  small  for  material  2 as  compared  to 
material  1.  However,  change  in  system  loss  factor  with  the 
application  of  damping  tape  for  both  systems  is  about  the 
same.  For  material  2,  neglecting  the  base  structure  loss 
factors  would  yield  reasonable  results  for  system  response, 
but  for  material  1 (as  is  true  for  most  composite  materials) 
the  base  structure  loss  factor  cannot  be  neglected. 

The  variation  of  tip  displacement  with  the  amount  of 
treatment  for  materials  1 and  2 is  shown  in  Figures  8.7 
through  8.10.  For  both  materials,  the  vibration  amplitude  is 
dramatically  reduced  with  increasing  amounts  of  damping 
treatment.  As  can  be  seen  from  the  damping  curves,  the 
presence  of  an  optimal  length  and  location  for  a given 
structure  and  loading  condition  is  evident.  A clear  optimum 
length  for  the  damping  treatment  is  seen  from  the  damping 
curves,  but  the  tip  displacement  variation  does  not  show  the 
same  clear  result.  This  is  because  displacement  is  function 
of  both  damping  and  stiffness.  It  can  also  be  seen  that 
treatment  closer  to  the  root  of  the  beam  will  have  a stronger 
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Figure  8.6.  Analytical  and  Experimental  Results  for 
Material  2 a/L  = 0.0. 
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Figure  8.7.  Variation  of  Tip  Displacement 
(material  1,  a/L  = 0.0) 
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Figure  8.8. 


Effect  of  Location  on  Tip  Displacement 
(material  1,  mode  1) 
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Figure  8.9. 


Effect  of  Location  on  Tip  Displacement 
(material  1,  mode  2) 
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Figure  8.10.  Variation  of  Tip  Displacement 
(material  2,  a/L  = 0.0) 
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influence  on  system  response  for  mode  1,  while  treatment 
located  away  from  the  root  of  the  beam  has  a stronger 
influence  on  mode  2 response.  Several  other  factors  such  as 
the  thickness  of  the  damping  layer  and  the  type  of  damping 
treatment  also  affect  the  system  response  strongly. 

8.2  Effect  of  Prestress  on  Loss  Factor  and  Frequency 

Structural  damping  of  a composite  beams  subjected  to  an 
axial  preload  was  evaluated  by  using  a finite  element 
approach,  and  some  of  the  results  were  compared  with  the 
experimental  results  of  Mantena  [87].  Results  of  the  effects 
of  different  parameters,  such  as  length  of  the  treatment  and 
structural  prestress,  on  the  overall  damping  of  the  composite 
system  are  presented.  The  specimen  used  was  a five  ply 
graphite  epoxy  beam  of  stacking  sequence  [90/90/0/90/90], 
length  0.2032  m,  width  0.01905  m,  and  thickness  of  0.00064  m. 
The  material  properties  of  the  three  layers  are  given  in  Table 
8.2. 

Figures  8.11  to  8.16  show  the  effect  of  preload  on  loss 
factor  and  frequency  for  different  lengths  of  the  damping 
treatment.  In  all  cases  the  damping  tape  is  applied  starting 
at  the  fixed  end  of  the  beam.  Frequency  dependent  properties 
of  the  damping  material  at  room  temperature  were  used  in  the 
calculations.  For  the  static  part  of  the  loading,  the  beam  is 
fixed  at  one  end  with  an  axial  load  acting  at  the  free  end, 
and  for  the  dynamic  part  of  the  loading,  the  beam  is  fixed  at 
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Table  8.2.  Material  Properties 


Graphite 

Epoxy 

VISCOELASTIC 

CORE 

CONSTRAINING 
LAYER  AND  A1 
FACING 

P (g/cm3) 

1.58 

0.980 

2 .76 

EL(GPa) 

127.90 

Fig.  2.3 

69 

Et  ( GPa ) 

10.27 

- 

69 

Glt  (GPa) 

7 .31 

- 

26 

VLT 

0.22 

0.499 

0.32 

‘Hl 

0.003 

Fig . 2.3 

0.005 

rij 

0.015 

Fig.  2.3 

0 . 005 
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Figure  8.11. 


Variation  of  Frequency  with  Preload 
(a/L  = 0,  b/L  = 0.2) 


LOSS  FACTOR 
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Figure  8.12.  Variation  of  Loss  Factor  with  Preload 
(a/L  = 0,  b/L  = 0.2) 


FREQUENCY  (Hz) 
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Figure  8.13. 


Variation  of  Frequency  with  Preload 
(a/L  = 0,  b/L  = 0.4) 


LOSS  FACTOR 
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Figure  8.14. 


Variation  of  Loss  Factor  with  Preload 
(a/L  = 0,  b/L  = 0.4) 
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Figure  8.15. 


Variation  of  Frequency  with  Preload 
(a/L  = 0,  b/L  = 0.6) 


LOSS  FACTOR 
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Figure  8.16.  Variation  of  Loss  Factor  with  Preload 
(a/L  = 0,  b/L  = 0.6) 
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both  ends.  In  all  the  cases  damping  decreases  with  increasing 
preload  for  each  of  the  first  three  modes  but  the  stiffness  of 
the  structure  has  increased  at  the  same  time.  Decreased  loss 
factor  alone,  therefore,  does  not  necessarily  result  in 
increased  vibration  amplitude  (Figure  8.17)  since  dynamic 
displacement  depends  on  both  stiffness  and  damping.  However, 
for  a given  preload,  the  displacement  is  always  considerably 
smaller  when  the  damping  treatment  is  present. 

Figures  8.18  to  8.21  show  the  results  of  loss  factor  and 
freguency  of  fully  taped  and  bare  beams  for  different  preloads 
compared  with  experimental  results  of  Mantena  et  al.  At  the 
preload  of  45  N,  the  analytical  and  experimental  first  mode 
freguencies  are  significantly  different.  At  other  preloads 
the  results  of  both  freguency  and  loss  factor  are  in  good 
agreement  before  any  failure  is  initiated.  Failure  is 
believed  to  be  initiated  at  the  point  at  which  the  loss  factor 
shows  an  increase  with  increasing  preload  in  the  experimental 
result.  At  high  preloads  the  difference  between  finite 
element  results  and  experimentally  measured  loss  factors  is 
due  to  the  problems  associated  with  measurement  of  small 
levels  of  damping. 
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Figure  8.17.  Response  Spectrum  for  Different  Tape  Lengths 
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Figure  8.18. 


Variation  of  Frequency  with  Preload 
(a/L  = 0,  b/L  = 1.0) 


LOSS  FACTOR 
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Figure  8.19.  Variation  of  Loss  Factor  with  Preload 
(a/L  = 0,  b/L  = 1.0) 
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Figure  8.20.  Effect  of  Compressive  Preload  on  Frequency 
(a/L  = 0,  b/L  = 1.0) 
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Figure  8.21.  Effect  of  Compressive  Preload  on  Loss  Factor 
(a/L  = 0,  b/L  = 1.0) 
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8 . 3 Rotating  Laminated  Tapered  Beam 

The  elements  developed  in  the  present  work  are  used  to 
study  the  effect  of  rotation  on  the  free  vibration  frequencies 
modal  loss  factors  and  mode  shapes  of  a laminated  tapered 
beam.  The  properties  and  dimensions  of  the  tapered  beam  are 
given  in  Figure  8.22. 

The  variation  of  frequencies  and  loss  factors  with 
rotation  speeds  for  both  the  bare  beam,  and  the  beam  fully 
treated  with  the  viscoelastic  damping  tape,  for  three 
different  rotation  speeds,  are  given  in  Tables  8.3  and  8.4. 
Frequency-invarient  properties  given  in  Table  7.2  were  used 
for  the  viscoelastic  material.  A comparison  of  the  vibration 
modes  shapes  of  the  fixed  beam,  and  the  beam  rotating  at  500 
rad/sec  is  made  in  Figures  8.23  to  8.30.  These  results  are 
for  the  untreated  beam.  Results  show  that  a large  shift  in 
resonant  frequencies  is  caused  by  the  centrifugal  stiffening 
effect.  The  vibration  modes  and  mode  shapes  are  also  at 
different  relative  positions  for  the  two  cases. 

Since  the  relative  strain  energy  distribution  can  be 
modified  considerably  by  prestress  it  is  important  to  consider 
prestress  effects  in  designing  passive  damping  treatments.  In 
the  present  case,  a damping  design  that  is  effective  for  the 
nonrotating  beam  is  not  necessarily  the  best  design  for  the 
rotating  beam. 
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Ply  thickness  0.000127  mm 
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Length 

1 m 

e2 

= 10.27  GPa 

Width 

0.05  m (fixed  end) 

Gia 

= 7.31  GPa 

0.02  m (free  end) 

1)12 

^23 
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Figure  8.22.  Rotating  Laminated  Tapered  Beam 
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Table  8.3.  Results  for  the  Rotating  Tapered  Beam 


MODE  NUMBER 

n 

(rad/sec) 

FREQUENCY 

(Hz) 

LOSS 

FACTOR 

1 

000 

4.81 

0.00500 

1 

100 

17.63 

0.00053 

1 

500 

81.55 

0 . 00008 

2 

000 

25.36 

0 . Q0500 

2 

100 

45.03 

0.00161 

2 

500 

121.10 

0.00231 

3 

000 

67 . 87 

0.00500 

3 

100 

85.46 

0.00476 

3 

500 

176.00 

0.00198 

4 

000 

81.17 

0.00500 

4 

100 

90.24 

0.00285 

4 

500 

184 . 51 

0.00014 

5 

000 

108.98 

0.00500 

5 

100 

112.53 

0.00469 

5 

500 

296.41 

0.00035 

6 

000 

132 . 58 

0.00500 

6 

100 

156.58 

0.00360 

6 

500 

355.20 

0.00285 

7 

000 

221.21 

0 . 00500 

7 

100 

246 . 02 

0.00405 

7 

500 

415.87 

0.00386 

8 

000 

262 . 12 

0 . 00500 

8 

100 

266.80 

0.00483 

8 

500 

422 . 83 

0.00064 
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Table  8.4.  Results  for  the  Damped  Rotating  Tapered  Beam 


MODE  NUMBER 

D (rad/sec) 

FREQUENCY 

(Hz) 

LOSS  FACTOR 

1 

000 

4 . 84 

0.00736 

1 

100 

17 . 65 

0.00678 

1 

500 

81.74 

0.00063 

2 

000 

25.47 

0.01786 

2 

100 

45.13 

0.02122 

2 

500 

116.84 

0.00239 

3 

000 

68 . 10 

0.02977 

3 

100 

77 . 60 

0.00502 

3 

500 

175.13 

0 . 04074 

4 

000 

75.47 

0.00527 

4 

100 

90.52 

0.04412 

4 

500 

184 . 99 

0.00115 

5 

000 

107 .38 

0.11155 

5 

100 

110.98 

0.14983 

5 

500 

297.12 

0.00271 

6 

000 

133.30 

0.03998 

6 

100 

157 .43 

0.06415 

6 

500 

350.36 

0 . 04747 

7 

000 

224 . 03 

0.04594 

7 

100 

249 . 00 

0.07841 

7 

500 

376.33 

0.04769 

8 

000 

255.40 

0.09314 

8 

100 

260.17 

0.13789 

8 

500 

424 . 87 

0 . 00548 
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Figure  8.23.  Vibration  Mode  1:  Laminated  Tapered  Beam 
(Top  f2  = 0,  Bottom  n = 500) 
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Figure  8.24.  Vibration  Mode  2:  Laminated  Tapered  Beam 
(Top  n = 0,  Bottom  n = 500) 
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Figure  8.25.  Vibration  Mode  3:  Laminated  Tapered  Beam 
(Top  n = 0,  Bottom  n = 500) 
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Figure  8.26.  Vibration  Mode  4:  Laminated  Tapered  Beam 
(Top  n = 0,  Bottom  n = 500) 


135 


Figure  8.27.  Vibration  Mode  5:  Laminated  Tapered  Beam 
(Top  ft  = 0,  Bottom  ft  = 500) 
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Figure  8.28.  Vibration  Mode  6:  Laminated  Tapered  Beam 
(Top  n = 0,  Bottom  n = 500) 
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Figure  8.29.  Vibration  Mode  7:  Laminated  Tapered  Beam 
(Top  n = 0,  Bottom  n = 500) 
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Figure  8.30. 


Vibration  Mode  8:  Laminated  Tapered  Beam 
(Top  ft  = 0,  Bottom  ft  = 500) 
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8.4  Laminated  Rotating  Circular  Annulus 

In  this  section,  the  effect  of  rotation  on  the  free 
vibration  frequencies  and  mode  shapes  of  an  internally 
clamped,  laminated  circular  annulus  are  presented.  The 
annulus  has  an  inner  radius  of  0.1  m and  an  outer  radius  of 
0.5  m.  The  material  properties  are  the  same  as  that  for  the 
example  in  Section  8.3,  and  the  stacking  sequence  is  [0/90]5s. 

The  first  few  resonant  frequencies  and  loss  factors  for 
the  disk  rotating  at  0 and  500  rpm  are  given  in  Table  8.5. 
Similar  results  for  the  disk  with  damping  treatment  applied 
fully  to  one  face  of  the  disk  are  given  in  Table  8.6.  The 
damping  layer  was  assumed  to  have  the  frequency-independent 
properties  given  in  Table  7.2.  The  effect  of  the  prestress  on 
the  mode  shapes  of  the  bare  annulus  is  shown  in  Figures  8.31 
to  8.36.  Six  eight-node  plate  elements  were  used  in  both  the 
radial  and  the  circumferential  directions  in  the  finite 
element  model. 

As  in  previous  cases,  the  stiffening  due  to  rotation 
causes  an  increase  in  frequency  and  a corresponding  reduction 
in  loss  factor.  The  drop  in  loss  factor  for  each  mode  is  not 
uniform  and  depends  on  how  much  the  strain  energy  distribution 
in  a mode  is  altered  due  to  the  prestress.  In  each  case,  the 
damping  ability  of  the  disk  is  considerably  higher  when  the 
damping  layer  is  present. 
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Table  8.5.  Results  for  the  Rotating  Circular  Annulus 


MODE  NUMBER 

n 

(rad/sec) 

FREQUENCY 

(Hz) 

LOSS 

FACTOR 

1 

000 

13 . 08 

0.00500 

1 

500 

75.48 

0.00034 

2 

000 

14 .30 

0.00500 

2 

500 

89 . 44 

0.00017 

3 

000 

22 .71 

0.00500 

3 

500 

121.14 

0 . 00018 

4 

000 

38 .93 

0.00500 

4 

500 

157 . 03 

0.00031 

5 

000 

70.45 

0.00500 

5 

500 

200.44 

0.00062 

6 

000 

93 .83 

0 . 00500 

6 

500 

234.62 

0.00084 
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Table  8.6.  Results  for  the  Damped  Rotating  Circular  Annulus 


MODE 

NUMBER 

n (rad/sec) 

FREQUENCY 

(Hz) 

LOSS  FACTOR 

1 

000 

15.09 

0.04421 

1 

500 

77 . 33 

0.00361 

2 

000 

16.20 

0.03090 

2 

500 

90.83 

0.00212 

3 

000 

26.95 

0.02500 

3 

500 

122 .83 

0.00111 

4 

000 

51.14 

0.03238 

4 

500 

162 . 06 

0.00307 

5 

000 

97 .30 

0.05525 

5 

500 

218 . 12 

0.00911 

6 

000 

106.32 

0.06906 

6 

500 

242 . 45 

0.01392 
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Figure 


.31.  Vibration  Mode  1:  Laminated  Circular  Annulus 
(Top  Cl  = 0,  Bottom  n = 500) 
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Figure 


.32.  Vibration  Mode  2:  Laminated  Circular  Annulus 
(Top  n = 0,  Bottom  n = 500) 
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Figure  8.33. 


Vibration  Mode  3:  Laminated  Circular  Annulus 
(Top  n = 0,  Bottom  f2  = 500) 
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Figure  8.34.  Vibration  Mode  4:  Laminated  Circular  Annulus 
(Top  n = 0,  Bottom  n = 500) 
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Figure  8.35. 


Vibration  Mode  5:  Laminated  Circular  Annulus 
(Top  ft  = 0,  Bottom  ft  = 500) 
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Figure  8.36. 


Vibration  Mode  6:  Laminated  Circular  Annulus 
(Top  n = 0,  Bottom  n = 500) 


CHAPTER  9 

SUMMARY,  CONCLUSIONS  AND  SUGGESTIONS  FOR  FUTURE  WORK 

Dynamic  response  characteristics  of  structures  can  be 
improved  significantly  by  capitalizing  on  the  high  damping 
ability  of  some  polymeric  materials.  It  is  known  that  by 
incorporating  the  high  damping  polymeric  materials  in  a 
constrained  form,  the  resonant  noise  and  vibration  problems  in 
complex  structures  can  be  suppressed  effectively.  However,  in 
order  to  optimize  the  controlled  energy  dissipation 
characteristics  of  structures,  accurate  analysis  capability  is 
required. 

A major  objective  of  this  research  was  to  examine  the 
changes  in  the  damping  properties  of  composite  structural 
elements  containing  viscoelastic  damping  layers.  Analytical 
capabilities  to  predict  the  total  damping  in  generally 
laminated  composites,  including  both  internal  damping  and  a 
component  due  to  the  damping  treatments,  were  developed. 
Specialized  elements,  and  techniques  for  the  estimation  of 
overall  structural  loss  factor,  a measure  of  damping,  were 
formulated  and  implemented  in  a stand-alone  finite  element 
code  UFPAC . The  program,  UFPAC,  was  developed  entirely  in 
this  study.  Particular  attention  was  focused  on  the 
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efficiency  of  the  formulations.  Analysis  capabilities  include 
modal  damping  analysis  using  the  modal  strain  energy  approach, 
and  forced  harmonic  vibration  analysis,  for  generally 
laminated  composite  structures  about  prestressed 
configurations . 

On  the  basis  of  the  results  presented,  it  can  be  seen 
that  the  finite  element  model  developed  in  this  study,  can 
accurately  and  efficiently  represent  generally  laminated, 
viscoelastically  damped  composites.  The  results  also  show 
that  factors  like  the  relative  stiffnesses  of  the  different 
layers,  location  and  amount  of  the  damping  layer,  and  the 
distribution  of  the  damping  layer  can  be  altered  to  construct 
integrally  damped,  optimally  configured,  laminated  composite 
structural  elements.  Initial  stress  is  shown  to  have  a 
pronounced  effect  on  the  vibration  characteristics  of  a 
structure;  therefore,  it  is  important  to  consider  the  effect 
of  initial  stress  in  designing  a passively  damped  structure. 
The  results  presented  in  Chapter  8 show  that  if  the  prestress 
has  a large  stiffening  effect,  it  may  be  difficult  to  improve 
vibration  characteristics  by  using  passive  damping  materials. 
This  is  because  the  elastic  stiffness  of  the  structure  is 
dominated  by  the  stiffening  effect  of  the  preload.  It  can  be 
seen  from  the  development  presented  in  Chapter  5 that  the 
increase  in  stiffness  due  to  preload  contributes  only  to 
energy  stored  and  not  the  energy  dissipated,  and  therefore, 
causes  a reduction  in  overall  damping  ability. 
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Additional  work  is  still  required  in  the  area  of  using 
the  analysis  methods  as  a practical  design  tool.  Since  a 
number  of  variables  have  to  be  considered  in  designing 
passively  damped  composite  structures,  the  resulting 
integrated  design  problem  is  quite  complicated.  The  work 
presented  in  this  study  has  been  used  in  the  analysis  phase  of 
some  optimization  studies  dealing  with  the  problem  of 
integrally  damped  composite  structural  elements  [87,88]. 

Some  applications  require  high  damping  ability  at  very 
low  frequencies  and  deformation  levels.  Experimental 
evidence,  and  modifications  in  analysis  methods  are  needed 
with  specific  reference  to  performance  of  the  viscoelastic 
damping  materials  in  such  situations. 


APPENDIX  A 

COMBINED  STIFFNESS  COEFFICIENTS  ABOUT  AN  ARBITRARY  PLANE 

The  formulation  of  the  laminated  beam  and  plate  elements 
presented  in  Chapter  4 required  the  combined  stiffnesses 
specified  about  an  arbitrary  reference  plane.  Some  sample 
calculations  of  laminate  combined  stiffnesses  about  different 
reference  planes,  are  presented  in  this  appendix. 

In  the  first  example,  a symmetrically  stacked  laminate  is 
considered,  and  in  the  second  example,  a generally  stacked 
laminate  is  considered.  The  following  calculations  are 
presented  for  each  case: 

(a)  Reference  plane  is  the  same  as  the  midplane 
(offset  = 0) 

(b)  Reference  plane  is  the  +z  surface  of  the  beam 
(offset  = -thickness/2) 

Graphite  epoxy  properties  from  Table  8.2  are  used.  For  the 
symmetric  stacking  sequence,  the  bending-stretching  coupling 
stiffness  coefficients  are  zero  if  the  mid-plane  is  specified 
as  the  reference  plane,  and  are  non-zero  if  the  plane  is 
offset.  In  the  second  example,  where  the  laminate  is 
generally  stacked,  the  stiffness  coefficient  matrix  is  fully 
populated  in  each  case.  The  stiffness  coefficient  submatrices 
are  in  the  order  [A],  [B],  [D]  and  [G] . 
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Example  1:  case  (a),  [0,15,30,45,60,75,90] 


103473868.2444691 

24332795.16948141 

28006103 . 11822166 


45.90135022686075 
5.499461208367602 
8 .213571478178362 

8582206.688065511 

1196625.671096968 


24332795.16948141 

103473870.8440072 

28006102.80990289 


5.499461208367602 

10.43858806521175 

4.316182211411334 

1196625.671096968 

8582206.743601156 


28006103 . 11822166 

28006102 .80990289 
32927197 .42326303 

1. 13686837721E-13 
-1.49213974509E-13 
-1. 27897692436E-13 

8.213571478178362 
4.316182211411334 
7 .763573402904583 


4 .547473508864E-13  -3 . 55271367880E-14 
3.552713678800E-14  -6 . 82121026329E-13 
1. 136868377216E-13  -1 . 49213974509E-13 


Example  1:  case  (b) , [0,15,30,45,60,75,90] 


103473868.2444691 

24332795.16948141 

28006103 . 11822166 

91988.26886933304 

21631.85490566897 

24897.42567209906 

127.6789212516978 

24.73018021950732 

30.34738290067442 

8582206.688065511 

1196625.671096968 


24332795.16948141 
103473870.8440072 

28006102 .80990289 

21631.85490566897 

91988.27118032237 

24897.42539800367 

24.73018021950732 

92.21616114451834 

26.44999339023660 

1196625.671096968 

8582206.743601156 


28006103 . 11822166 

28006102 .80990289 
32927197 . 42326303 

24897.42567209906 
24897.42539800367 
29272 . 27850928083 

30.34738290067442 

26.44999339023660 

33.78662899765524 
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Example  2:  case  (a),  [ 0 , 15 , 30 , 45 , 60 , 75 , 90 ] 2 


103473868.2444691  24332795.16948141 
24332795.16948141  103473870.8440072 
28006103.11822166  28006102 .80990289 


-19945.31077743465 
-1.01731515222E-04 
-2192 . 007599810531 


-1. 01731515222E-04 
19945.31098089768 
2192 . 007136977530 


28006103 . 11822166 

28006102 . 80990289 
32927197.42326303 

-2192 . 007599810531 
2192.007136977530 
-1. 01731515258E-04 

6.264876721946800 

6.264876556184358 

7.763573312465266 


28.16996894572135  5.499461117928285 
5.499461117928285  28.16996952722979 
6.264876721946800  6.264876556184358 

8582206.688065511  1196625.671096968 
1196625.671096968  8582206.743601156 


Example  2:  case  (b) , [ 0 , 15 , 30 , 45 , 60 , 75 , 90 ] 2 


103473868.2444691 
24332795.16948141 

28006103 . 11822166 

72042.95809189839 

21631.85480393746 

22705.41807228853 

74.48477740827961 

24.73017994818937 

24.50129863197974 

8582206.688065511 

1196625.671096968 


24332795.16948141 

103473870.8440072 

28006102.80990289 

21631.85480393746 
111933 . 5821612201 
27089.43253498120 

24.73017994818937 

145.4103055305725 

32.29607642455567 

1196625.671096968 

8582206.743601156 


28006103 . 11822166 

28006102 .80990289 
32927197.42326303 

22705.41807228853 
27089 .43253498120 
29272 .27840754932 

24.50129863197974 

32.29607642455567 

33.78662872633729 


APPENDIX  B 

SUBROUTINES  FOR  CALCULATING  ELEMENT  MATRICES  OF  THE  LAMINATED 

PLATE  ELEMENT 

In  this  section,  the  Fortran  subroutines  for  calculating 
the  finite  element  matrices  of  the  laminated  plate  element  are 
presented.  The  material  property  information  required  for  the 
calculation  is  the  combined  stiffness  matrix  described  in 
Appendix  A.  This  information  is  read  from  the  pre-prepared 
data  files  ABD1.DAT,  ABD2.DAT  etc.,  as  specified  in  the  input 
command  file  (Appendix  D) . The  other  aspects  of  the 
subroutines,  and  the  function  of  each  module  are  described  in 
the  subroutine  listing. 

The  subroutines  are  listed  in  the  following  order, 

STIFF6 (NN) : Calculates  the  element  stiffness  matrix 
of  element  number  NN. 

STDM6 (NN) : Calculates  the  values  of  shape  function, 
derivatives  of  shape  function,  Jacobian,  and  the 
strain  displacement  transformation  matrix  at  the 
Gauss  points.  The  subroutine  is  called  by  each  of 
the  other  subroutines. 

MASS6(NN):  Calculates  the  element  mass  matrix  of 

element  number  NN. 

GEOM_ST06 (NN) : Calculates  the  geometric  stiffness 

matrix  of  element  number  NN. 
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c 

c********************************************************** 

c 

C SUBROUTINE  EVALUATES  THE  ELEMENT  STIFFNESS  MATRIX 
C OF  ELMT  (A  VARIABLE  NODE  LAMINATED  OFFSET 
C PLATE  ELEMENT)  IN  THE  ELEMENT  COORDINATE  SYSTEM 
C 

C* ********************************************************* 

C 

SUBROUTINE  STIFF6 (NN) 

IMPLICIT  REAL  *8  (A-H,0-Z) 

PARAMETER  (IXX=03000 , IYY=0250) 

CHARACTER  FILE*5 

COMMON/EINFO/NTOT , NETY (0750), NPR (0750) 

COMMON/DATPROC/ID ( 6 , 01500) , ICONN ( 2 0 , 0750 ) , NEQ 
COMMON/ NODE/ NT, COOR( 01500, 3) 

COMMON/ESTIF/ES (120, 120) 

COMMON/E LMT/H (20) , DH ( 20 , 3 ) , B ( 10 , 12 0 ) , DET , NMAX , ETA , ZI , SI 
1 , A J (3,3) , AJI (3,3) 

COMMON/ WORKAREA/ ABD ( 8 , 8 ) , ES I ( 8 ) 

COMMON/ GAUSS/XG (4,4), WGT (4,4) 

C 

IF (NPR (NN) • LT. 10) THEN 

FILE=  ' ABD'//CHAR(48+NPR(NN) ) 

ELSE 

FILE=' ABD' //CHAR (48+NPR(NN) /10) //CHAR (48+NPR(NN) 

1- (NPR(NN) /10) *10) 

ENDIF 

OPEN  (UNIT=7 , NAME=FILE , STATUS= ' OLD ' ) 

REWIND  (UNIT  = 7) 

READ  IN  ABD  MATRIX 

DO  1=1,8 
DO  J=1 , 8 
ABD ( J , I ) =0 . DO 
ENDDO 
ENDDO 
C 

DO  I = 1,3 
DO  J=  1,3 

IF  (I.EQ.l)  THEN 

READ ( 7 , * ) (ABD ( J , K) , K=1 , 3 ) 

ELSE  IF  (I.EQ.2)  THEN 

READ ( 7 , * ) (ABD ( J , K) , K=4 , 6 ) 

DO  K =1,3 

ABD ( J+3 , K)  = ABD ( J , K+3 ) 

ENDDO 

ELSE 

READ (7, *) (ABD (J+3 , K) ,K=4, 6) 

ENDIF 

ENDDO 
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ENDDO 

READ ( 7 , * ) ABD (7,7) , ABD (7,8) 

READ ( 7 , * ) ABD (8,7), ABD (8,8) 

C 

C CHECK  TO  SEE  THE  NODES  PRESENT  IN  THE  VARIABLE 
C NODE  ELEMENT 
C 

DO  1=1,9 

IF ( ICONN ( I , NN) . NE . 0) THEN 
NMAX=I 
ENDIF 

ENDDO 

C 

DO  1=1 , NMAX*6 
DO  J=1 , NMAX*6 
ES ( J, I) =0 . DO 
ENDDO 
ENDDO 
C 

C CALCULATE  THE  ELEMENT  STIFFNESS  MATRIX 

C NINT  ORDER  OF  INTEGRATION 

C 

NINT=NETY (NN) - (NETY (NN) /10) *10 
DO  LX=1 , NINT 

ETA=XG (LX, NINT) 

C 

DO  LY=1 , NINT 
ZI=XG ( LY , NINT) 

EVALUATE  DET  OF  JACOBIAN  AND  MATRIX  B AT  (ETA,ZI) 
CALL  STDM6 (NN) 

ADD  CONTRIBUTION  TO  THE  ELEMENT  STIFFNESS  MATRIX 

WT=WGT ( LX , NINT ) *WGT ( LY , NINT) *DET 
C 

DO  1=1 , NMAX*6 
DO  K=1 , 8 
ESI (K) =0 
DO  J=1 , 8 

ESI (K) =ESI (K)+B(J, I) *ABD(J,K) 

ENDDO 

ENDDO 

C 

DO  L=I,NMAX*6 
STIFF=0 . DO 
DO  M=1 , 8 

STIFF=STIFF+ESI (M) *B (M, L) 

ENDDO 

ES ( I , L) =ES ( I , L) +STIFF*WT 
ENDDO 
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C 


C 


C 


ENDDO 


ENDDO 

ENDDO 

DO  1=1 , NMAX*6 

DO  J=I+1 , NMAX*6 
ES (J, I) =ES (I , J) 
ENDDO 
ENDDO 

RETURN 

END 
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C 

c******************************************************* 

c 

C SUBROUTINE  TO  EVALUATE  STRAIN-DISPLACEMENT 
C TRANSFORMATION  AT  A POINT  (ETA,ZI)  OF  AN 
C ISOPARAMETRIC  VARIABLE  NODE  PLATE  ELEMENT  (ELMT) 

C 

c 

SUBROUTINE  STDM6 (NN) 

IMPLICIT  REAL  *8  (A-H,0-Z) 

PARAMETER  ( IXX=03 000 , IYY=0250 ) 

COMMON/DATPROC/ID (6,01500) , ICONN (20 , 0750) , NEQ 
COMMON/NODE/NT, COOR( 015 00, 3) 

COMMON/ ELMT/H (20) , DH (20 , 3 ) , B ( 10 , 120) , DET , NMAX , ETA , ZI , SI 
1 , A J (3,3) , AJI (3,3) 

INTERPOLATION  FUNCTIONS 

H ( 1 ) = 0.25D0*(1+ZI)*( 1-ETA) 

H ( 2 ) = 0 . 25D0* ( 1-ZI) * ( 1-ETA) 

H ( 3 ) = 0.25D0* (1-ZI) * (1+ETA) 

H ( 4 ) = 0.25D0*(1+ZI)*( 1+ETA) 

IF ( ICONN ( 5 , NN) . EQ. 0) THEN 
H ( 5 ) = 0 . DO 
ELSE 

H ( 5 ) = 0.5D0*(1-ZI**2) * ( 1-ETA) 

DO  1=1,2 

H ( I ) = H ( I ) - 0 . 5D0*H ( 5 ) 

ENDDO 
ENDIF 

IF (ICONN ( 6 , NN) . EQ. 0) THEN 
H ( 6 ) = 0 . DO 
ELSE 

H ( 6 ) = 0. 5D0* (1-ETA**2) * (1-ZI) 

DO  1=2,3 

H ( I ) = H ( I ) - 0 . 5D0*H ( 6 ) 

ENDDO 
ENDIF 

I F ( I CONN ( 7 , NN ) . EQ. 0) THEN 
H ( 7 ) = 0 . DO 
ELSE 

H ( 7 ) = 0. 5D0* (1-ZI**2) * (1+ETA) 

DO  1=3,4 

H ( I ) = H ( I ) - 0 . 5D0*H ( 7 ) 

ENDDO 
ENDIF 

I F ( I CONN ( 8 , NN ) . EQ. 0) THEN 
H ( 8 ) = 0 . DO 
ELSE 

H ( 8 ) = 0 . 5 DO* ( 1-ETA**2 ) * ( 1+ZI ) 

DO  1=1, 4, 3 
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H ( I ) = H ( I ) - 0 . 5D0*H ( 8 ) 

ENDDO 

ENDIF 

IF ( ICONN ( 9 , NN) . EQ . 0 ) THEN 

H ( 9 ) = 0 . DO 
ELSE 

H ( 9 ) = ( 1 . D0-ETA*ETA) * (1 . DO-ZI*ZI) 

DO  1=1,4 

H ( I ) = H ( I ) + 0 . 25D0*H (9) 

H ( 1+4 ) = H ( 1+4 ) - 0 . 5D0*H ( 9 ) 

ENDDO 

ENDIF 

C 

C DERIVATIVES  OF  INTERPOLATION  FUNCTIONS  W.R.T. 

C THE  NATURAL  COORDINATES 

C 

C WITH  RESPECT  TO  ETA  (LOCAL  X COORDINATE) 

C DERIVATIVE  OF  H WRT  ETA 

C 

DH (1,1)=  -0 • 25D0* ( 1+ZI ) 

DH (2,1)=  -0. 25D0* (1-ZI) 

DH (3,1)=  0. 25D0* (1-ZI) 

DH (4,1)=  0. 25D0* (1+ZI) 

I F ( I CONN ( 5 , NN ) . EQ. 0) THEN 

DH (5,1)=  0 . DO 
ELSE 

DH (5,1)=  (-0.5D0) *(1-ZI**2) 

DO  1=1,2 

DH ( I , 1 ) = DH (1,1)  - 0. 5D0*DH(5, 1) 

ENDDO 

ENDIF 

IF ( ICONN ( 6 , NN) . EQ. 0) THEN 

DH (6,1)=  0 . DO 
ELSE 

DH (6,1)=  -ETA* ( 1-ZI ) 

DO  1=2,3 

DH ( I , 1) = DH (1,1)  - 0. 5D0*DH(6, 1) 

ENDDO 

ENDIF 

IF (ICONN ( 7 , NN) . EQ. 0) THEN 

DH (7,1)=  0 . DO 
ELSE 

DH (7,1)=  0.5D0*(1-ZI**2) 

DO  1=3,4 

DH ( I , 1 ) = DH (1,1)  - 0 . 5D0*DH (7,1) 

ENDDO 

ENDIF 

I F ( I CONN ( 8 , NN ) . EQ. 0) THEN 

DH (8,1)=  0 . DO 
ELSE 

DH (8,1)=  -ETA* ( 1+ZI) 

DO  1=1, 4, 3 
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DH ( I , 1 ) = DH (1,1)  - 0 . 5D0*DH (8,1) 
ENDDO 
ENDIF 

IF ( ICONN ( 9 , NN) . EQ . 0 ) THEN 

DH (9,1)=  0 . DO 
ELSE 

DH (9,1)=  -2 . DO* ETA* (1.D0— ZI*ZI) 

DO  1=1,4 

DH ( I , 1 ) = DH (1,1)  + 0 . 25D0*DH (9,1) 

DH ( 1+4 , 1 ) = DH ( 1+4 , 1 ) - 0 . 5D0*DH (9,1) 
ENDDO 
ENDIF 
C 

C WITH  RESPECT  TO  ZI  (LOCAL  Y COORDINATE) 

C DERIVATIVE  OF  H WRT  ZI 
C 

DH (1,2)=  0.25D0* (1-ETA) 

DH (2,2)=  -0. 25D0* (1-ETA) 

DH (3,2)=  -0 . 25D0* ( 1+ETA) 

DH (4,2)=  0 . 2 5D0* ( 1+ETA) 

IF ( ICONN ( 5 , NN) • EQ. 0) THEN 

DH (5,2)=  0 . DO 
ELSE 

DH (5,2)=  -ZI* ( 1-ETA) 

DO  1=1,2 

DH ( I , 2 ) = DH (1,2)  - 0. 5D0*DH(5, 2) 
ENDDO 
ENDIF 

I F ( I CONN ( 6 , NN ) . EQ. 0) THEN 

DH (6,2)=  0 . DO 
ELSE 

DH (6,2)=  (-0.5D0) * (1-ETA**2) 

DO  1=2,3 

DH ( I , 2 ) = DH (1,2)  - 0 . 5D0*DH (6,2) 
ENDDO 
ENDIF 

I F ( I CONN ( 7 , NN ) . EQ. 0) THEN 

DH (7,2)=  0 . DO 
ELSE 

DH (7,2)=  -ZI* ( 1+ETA) 

DO  1=3,4 

DH ( I , 2 ) = DH (1,2)  - 0 . 5D0*DH (7,2) 
ENDDO 
ENDIF 

IF ( ICONN ( 8 , NN) . EQ . 0) THEN 

DH (8,2)=  0 . DO 
ELSE 

DH (8,2)=  0 . 5D0* ( 1-ETA**2 ) 

DO  1=1,4, 3 

DH ( I , 2 ) = DH (1,2)  - 0 . 5D0*DH (8,2) 
ENDDO 
ENDIF 
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IF ( ICONN ( 9 , NN) . EQ. 0) THEN 
DH (9,2)=  0 . DO 
ELSE 

DH (9,2)=  -2 . D0*ZI* ( 1 . D0-ETA*ETA) 

DO  1=1,4 

DH ( 1 , 2 ) = DH (1,2)  + 0 . 25D0*DH (9,2) 

DH ( 1+4 , 2 ) = DH ( 1+4 , 2 ) - 0 . 5D0*DH ( 9 , 2 ) 

ENDDO 
ENDIF 

EVALUATE  THE  JACOBIAN  AT  (ETA,ZI) 

DO  1=1,2 

BUF1=0 . DO 
BUF2=0 . DO 
DO  K=  1,9 

BUFl=BUFl+COOR ( ICONN (K, NN) , I) *DH(K, 1) 
BUF2=BUF2+COOR ( ICONN (K, NN) ,1) *DH(K,2) 
ENDDO 

A J ( 1 , I ) = BUF1 
A J ( 2 , I ) = BUF2 

ENDDO 

CALCULATE  DETERMINANT  OF  THE  JACOBIAN  MATRIX 
AT  POINT  ( ETA , Z I ) 

DET  = AJ  (2,2)  *AJ  (1,1)  -AJ  (1,2)  * AJ  (2,1) 

IF  (DET.LE.O)  THEN 
PRINT*,'  ELEMENT  NO  ' , NN , 

1 ' : DETERMINANT  HAS  NEGATIVE  VALUE ' 

ENDIF 

CALCULATE  INVERSE  OF  THE  JACOBIAN 

BUF1=  1.D0/DET 
AJI (1,1) =AJ (2,2) *BUF1 
AJI (2,2) =AJ (1,1) *BUF1 
AJI ( 1 , 2 ) =-AJ (1,2) *BUF1 
AJI (2,1) =-AJ (2,1) *BUF1 
C 

C CALCULATE  GLOBAL  DERIVATIVE  OPERATOR 

C 

C 

DO  1=1 , NMAX*6 
DO  J=1 , 8 
B(J, I) =0 . DO 
ENDDO 
ENDDO 
C 

DO  1=1 , NMAX 
C 

DO  J=1 , 2 
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C 


C 


C 


B(l, 6*1-5) = 
B(2, 6*1-4)= 
B(3, 6*1-5)= 
B(3, 6*1-4)= 


B(l, 6*1-5) 
B ( 2 , 6*1-4) 
B(3, 6*1-5) 
B ( 3 , 6*1-4) 


+ AJI (1, J) *DH ( I , J) 
+ AJI (2, J) *DH ( I , J ) 
+ AJI (2, J) *DH (I , J) 
4-  AJI  (1,  J)  *DH  ( I , J) 


ENDDO 

B(4, 6*1-2)= 
B(5, 6*1-1)= 
B(6, 6*1-2)= 
B(6, 6*1-1)= 
B(7, 6*1-3)= 
B(7, 6*1-2)= 
B(8, 6*1-3)= 
B(8, 6*1-1)= 


B(l, 6*1-5) 
B(2, 6*1-4) 
B ( 3 , 6*1-5) 
B(3, 6*1-4) 
B(l, 6*1-5) 
-H  ( I ) 

B(2, 6*1-4) 
“H  ( I ) 


ENDDO 


END  CALCULATION  OF  STRAIN  DISPLACEMENT  MATRIX 


RETURN 

END 
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C 

C********************************************************** 

c 

C SUBROUTINE  EVALUATES  THE  ELEMENT  MASS  MATRIX 

C OF  ELMT  (A  VARIABLE  NODE  LAMINATED  OFFSET 

C PLATE  ELEMENT)  IN  THE  ELEMENT  COORDINATE  SYSTEM 
C 

C* ********************************************************* 

C 

SUBROUTINE  MASS6(NN) 

IMPLICIT  REAL  *8  (A-H,0-Z) 

PARAMETER  (IXX=03000 , IYY=0250) 

CHARACTER  FILE *5 

COMMON/EINFO/NTOT , NETY ( 0750 ) , NPR ( 0750 ) 

COMMON/ DATPROC/ ID ( 6 , 01500 ) , I CONN ( 2 0 , 07 50 ) , NEQ 
COMMON/MATPRO/AMP1 (6,30) 

COMMON/NODE/NT, COOR ( 015 00, 3) 

COMMON/ESTIF/EM (120,120) , EKI ( 120 , 120) 

COMMON/ELMT/H (20) , DH ( 2 0 , 3 ) , B ( 10 , 120 ) , DET , NMAX , ETA, ZI , SI 
1 , A J (3,3) , AJI (3,3) 

COMMON/ WORKAREA/ ABD ( 8 , 8 ) , ES I ( 8 ) 

COMMON/ GAUSS/XG (4,4), WGT (4,4) 

C 

C CHECK  TO  SEE  THE  NODES  PRESENT  IN  THE  VARIABLE 
C NODE  ELEMENT 
C 

DO  1=1,9 

IF ( ICONN ( I , NN) . NE . 0 ) THEN 
NMAX=I 
ENDIF 
DO  J=1 , 6 
ABD ( J , I ) =0 . DO 
ENDDO 

ENDDO 

C 

C 

C AMP1 ( 3 , NPR (NN) ) DENSITY 

C AMP1 ( 4 , NPR (NN) ) THICKNESS  OF  THE  ELEMENT 

C AMP1 ( 5 , NPR (NN) ) NODE  OFFSET 

C 

RH0=AMP1 ( 1 , NPR (NN) ) 

TH=AMP1 ( 2 , NPR ( NN ) ) 

0FF=AMP1 ( 3 , NPR ( NN ) ) 

C 

ABD (1,1)=  RHO*TH 
ABD (1,4)=  -OFF*TH*RHO 
ABD (4,1)=  ABD (1,4) 

ABD (2,2)=  RHO*TH 
ABD (2,5)=  -OFF*TH*RHO 
ABD (5,2)=  ABD (2, 5) 

ABD (3,3)=  RHO*TH 

ABD (4,4)=  RHO* (OFF**2*TH+TH**3/12 . DO) 
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ABD (5,5)=  ABD (4,4) 

DO  1=1 , NMAX*6 
DO  J=1,NMAX*6 
EM ( J , I ) =0 . DO 
ENDDO 
ENDDO 

CALCULATE  THE  ELEMENT  STIFFNESS  MATRIX 
NINT  ORDER  OF  INTEGRATION 

NINT=NETY (NN) - (NETY (NN) /10) *10 
NINT=3 

DO  LX=1 , NINT 

ETA=XG ( LX , NINT ) 

DO  LY=1 , NINT 
ZI=XG ( LY , NINT) 

EVALUATE  DET  OF  JACOBIAN  AND  MATRIX  B AT  (ETA,ZI) 

CALL  STDM6 (NN) 

DO  1=1 , NMAX 
DO  J=1 , 6 
DO  K=1 , 6 
IF (K. EQ. J) THEN 
B ( K, (1-1) *6+J) =H ( I ) 

ELSE 

B (K, (1-1) *6+J) =0 . DO 
ENDIF 
ENDDO 
ENDDO 
ENDDO 

ADD  CONTRIBUTION  TO  THE  ELEMENT  STIFFNESS  MATRIX 

WT=WGT ( LX , N I NT ) *WGT ( LY , NINT) *DET 

DO  1=1 , NMAX*6 
DO  K=1 , 6 
ESI ( K) =0 . DO 
DO  J=1 , 6 

ESI (K) =ESI (K) +B ( J , I) *ABD ( J , K) 

ENDDO 

ENDDO 

DO  L=I , NMAX*6 
EMASS=0 . DO 
DO  M=1 , 6 

EMASS=EMASS+ESI (M) *B(M, L) 
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ENDDO 

EM ( I , L) =EM ( I , L) +EMASS*WT 
ENDDO 
ENDDO 

ENDDO 

ENDDO 

DO  1=1 , NMAX*6 

DO  J=I+1 , NMAX*6 
EM ( J , I ) =EM ( I , J ) 

ENDDO 

ENDDO 

RETURN 

END 
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C 

c* **************************************************** 
c 

C SUBROUTINE  CALCULATES  THE  GEOMETRIC  MATRIX 
C BY  USING  RESULTS  OF  A STATIC  ANALYSIS 
C 

C* **************************************************** 

c 

SUBROUTINE  GEOM_ST06 (NN) 

IMPLICIT  REAL  *8  (A-H,0-Z) 

PARAMETER (IXX=03 000, I YY=02 50) 

CHARACTER  FILE*5 

COMMON/E INFO/NTOT , NETY ( 0750 ) , NPR ( 0750 ) 

COMMON/ DATPROC/ ID (6,01500) , I CONN (20 , 0750) , NEQ 
COMMON/NODE/NT , COOR (01500,3) 

COMMON/ESTIF/ES (120,120) ,GF(IXX) 

COMMON/ELMT/H ( 20 ) ,DH(20,3) ,B(10,120) , DET, NMAX, ETA, ZI , SI 
1 , A J (3,3) , AJI (3,3) 

COMMON/ GAUSS/XG (4,4), WGT (4,4) 

COMMON/WORKAREA/A (3,3),BB(3,3),D(3,3),Z(3,3),Q0(2,2),Q1(2,2) 
1,Q2(2,2) , Q3 ( 2 , 2 ) ,E(8) ,ABD(9,9) ,S(18,18) , ESI (18) ,AN(6) 

1 , B1 (18,60) 

CALCULATE  GEOMETRIC  STIFFNESS  (INITIAL  STRESS 
STIFFNESS)  MATRIX 

IF (NPR (NN) .LT. 10) THEN 
FILE=  1 ABD ' //CHAR ( 48+NPR (NN) ) 

ELSE 

FILE- ' ABD ' //CHAR ( 4 8+NPR ( NN ) / 1 0 ) // CHAR ( 4 8 +NPR ( NN ) 

1- (NPR(NN) /10) *10) 

ENDIF 

OPEN  (UNIT=7 , NAME-FILE , STATUS- ' OLD ' ) 

REWIND  (UNIT  = 7) 

READ  IN  ABD  MATRIX 

DO  1=1,3 
DO  J=1 , 3 

IF  (I.EQ.l)  THEN 
READ ( 7 , * ) ( A ( J , K) , K=1 , 3 ) 

ELSEIF  (I.EQ.2)  THEN 
READ (7,*) ( BB ( J , K) ,K=1,3) 

ELSE 

READ (7,*) ( D ( J , K) ,K=1,3) 

ENDIF 
ENDDO 
ENDDO 
C 

READ ( 7 , * ) QO (1,1) , QO (1,2) 

READ ( 7 , * ) QO (2,1) ,Q0(2,2) 
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READ ( 7 , * ) Q1 (1,1) ,Q1(1,2) 

READ ( 7 , * ) Q1 ( 2 , 1 ) ,Q1(2,2) 

READ ( 7 , * ) Q2 (1,1) ,Q2(1,2) 

READ ( 7 , * ) Q2 (2,1) ,Q2(2,2) 

DO  1=1,3 

READ (7,*) (Z(I,J) , J=1 , 3 ) 

ENDDO 

READ ( 7 , * ) Q3 ( 1 , 1 ) ,Q3(1,2) 

READ ( 7 , * ) Q3 ( 2 , 1 ) ,Q3(2,2) 

C 

DO  1=1,8 
DO  J=l, 6 
ABD ( J , I ) =0  . DO 
ENDDO 
ENDDO 
C 

DO  1=1,9 

IF ( ICONN ( I , NN) . NE.O)  THEN 
NMAX=I 
ENDIF 
ENDDO 

INITIALIZE  ELEMENT  GEOMETRIC  STIFFNESS  MATRIX 

DO  1=1 , NMAX*6 
DO  J=1 , NMAX*6 
ES ( J , I) =0 . DO 
ENDDO 
ENDDO 

NINT=3 

START  CALCULATION  OF  GEOMETRIC  STIFFNESS  MATRIX 

DO  LX=1 , NINT 
ETA=XG ( LX , NINT ) 

C 

DO  LY=1 , NINT 
ZI=XG (LY , NINT) 

C 

CALL  STDM6 (NN) 

C 

DO  1=1 , NMAX*6 
DO  J=1 , 14 
B1 ( J, I) =0 . DO 
ENDDO 
ENDDO 
C 

DO  1=1 , NMAX 
C 


C 


DO  J=1 , 2 
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Bl(l, 6*1-5)  = B1  ( 1 , 6*1-5 ) + AJI (1, J) *DH(I, J) 
Bl(2, 6*1-5)  = B1 ( 2 , 6*1-5 ) + AJI (2 , J) *DH (I , J) 
B1 (4 , 6*1-4 ) = Bl(4, 6*1-4)  + AJI ( 1 , J ) *DH ( I , J) 
Bl(5, 6*1-4)  = Bl(5, 6*1-4)  + AJI ( 2 , J ) *DH ( I , J) 
B1 (7 , 6*1-3 ) = Bl(7, 6*1-3)  + AJI ( 1 , J ) *DH ( I , J ) 
B1  ( 8 , 6*1-3 ) = B1 ( 8 , 6*1-3 ) + AJI (2 , J) *DH (I , J) 
C 

Bl(10, 6*1-2)  = Bl(10, 6*1-2)  + AJI (1, J) *DH (I , J) 
Bl(ll, 6*1-2)  = B1 (11, 6*1-2)  + AJI (2 , J) *DH(I, J) 
Bl(13, 6*1-1)  =61(13,6*1-1)  + AJI  (1,  J)  *DH(I,  J) 
Bl(14, 6*1-1)  = Bl(14, 6*1-1)  +AJI(2,J)*DH(I,J) 
ENDDO 
C 

Bl(3, 6*1-2)  = -H ( I ) 

Bl(6, 6*1-1)  = -H ( I ) 

C 

ENDDO 

C 

WT=WGT ( LX , NINT ) *WGT ( LY , NINT) *DET 
C 

DO  1=1,8 
C 

E ( I ) =0 . DO 

C 

DO  J= 1 , NMAX 
C 

DO  K=1 , 6 

IF ( ID (K, I CONN ( J , NN) ) . NE . 0 ) THEN 
E ( I ) =E ( I ) +B ( I , (J-l) *6+K) *GF ( ID (K, ICONN ( J , NN) ) ) 
ENDIF 
ENDDO 
C 

ENDDO 

C 

ENDDO 

CALCULATE  KG=K1+K2+K3+K4 


DO  11=1,3 

CALCULATE  SUBMATRIX  K1 

IF  (II.EQ.l)  THEN 
DO  1=1,3 
DO  J=1 , 3 
ABD ( J , I ) =A ( J , I ) 

ABD ( J , 1+3 ) =BB ( J , I ) 

I F ( I . LT . 3 . AND . J . LT . 3 ) THEN 
ABD ( J+3 , 1+6 ) = Q0 ( I , J ) 
ENDIF 
ENDDO 
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ENDDO 

DO  1=1,5 
AN  (I)  =0  . DO 
DO  J=1 , 8 

AN ( I ) =AN ( I ) +ABD ( I , J ) *E(J) 

ENDDO 

ENDDO 

STRESS  RESULTANT  TENSOR  AT  THE  GAUSS  POINT  REPEATED  THREE 

ALONG  THE  DIAGONAL 

DO  1=1,9 

DO  J=l, 9 
S(I, J)=0.D0 
ENDDO 
ENDDO 

S ( 1 , 1 ) =AN ( 1 ) 

S ( 2 , 2 ) =AN ( 2 ) 

S ( 1 , 2 ) =AN ( 3 ) 

S ( 1 , 3 ) =AN ( 4 ) 

S ( 2 , 3 ) =AN ( 5 ) 

S (2 , 1) =S (1,2) 

S ( 3 , 2 ) =S (2,3) 

S ( 3 , 1) =S (1,3) 

DO  1=1,3 
DO  J=1 , 3 

S (J+3 , 1+3) =S (J, I) 

S ( J+6 , 1+6) =S ( J , I) 

ENDDO 

ENDDO 

CALCULATE  COMPONENT  DUE  TO  K2 

ELSE  IF ( II. EQ. 2) THEN 

DO  1=1,3 
DO  J=l, 3 

ABD ( J , I ) =-BB ( J , I ) 

ABD (J,I+3)=-D(J,I) 

IF ( I . LT . 3 . AND . J . LT . 3 ) THEN 
ABD (J+3 , 1+6 ) = Q1 ( I , J) 

ENDIF 

ENDDO 

ENDDO 

DO  1=1,5 

AN ( I ) =0 . DO 
DO  J=1 , 8 


C 


170 


AN ( I ) =AN ( I ) +ABD ( I , J ) *E(J) 

ENDDO 

ENDDO 

C 

C STRESS  RESULTANT  TENSOR  AT  THE  GAUSS  POINT  REPEATED  THREE 
TIMES 

C ALONG  THE  DIAGONAL 
C 

DO  1=10,18 
DO  J=1 , 9 
S ( J, I) =0 . DO 
ENDDO 
ENDDO 
C 

S ( 1 , 1 0 ) =-AN ( 1 ) 

S ( 2 , 11) =-AN ( 2 ) 

S ( 1 , 11 ) =-AN ( 3 ) 

S(l, 12)=-AN(4) 

S (2 , 12 ) =-AN (5) 

S (2 , 10) =S (1,11) 

S (3 , 11) =S (2,12) 

S ( 3 , 10) =S (1,12) 

C 

DO  1=10,12 
DO  J=l, 3 

S ( J+3 , 1+3 ) =S ( J , I) 

S ( J+6 , 1+3 ) =S ( J , I ) 

ENDDO 

ENDDO 

C 

C CALCULATE  COMPONENT  DUE  TO  K3 
C 

DO  1=10,18 
DO  J=1 , 9 

s (I , J)  =s  ( J,  I) 

ENDDO 

ENDDO 

C 

ELSE 

C 

DO  1=1,3 
DO  J=l, 3 
ABD(J, I) =D(J, I) 

ABD ( J , 1+3 ) =-Z (J,  I) 

IF ( I . LT . 3 . AND . J . LT . 3 ) THEN 
ABD (J+3,I+6)=  Q2 ( I , J) 

ENDIF 

ENDDO 

ENDDO 

C 

DO  1=1,5 
AN (I) =0. DO 


0 0 1-300 
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IMES 


DO  J=l,8 

AN ( I ) =AN ( I ) +ABD ( I , J ) *E(J) 

ENDDO 

ENDDO 

STRESS  RESULTANT  TENSOR  AT  THE  GAUSS  POINT  REPEATED  THREE 

ALONG  THE  DIAGONAL 

DO  1=9,18 

DO  J=9 , 18 
S ( I , J) =0 . DO 
ENDDO 
ENDDO 

S ( 1 0 , 1 0 ) =AN ( 1 ) 

S(ll, 11)=AN(2) 

S ( 1 0 , 1 1 ) =AN ( 3 ) 

S ( 10 , 12 ) =AN (4 ) 

S ( 11 , 12 ) =AN ( 5 ) 

S(11,10)=S(10,11) 

S ( 12 , 11) =S (11,12) 

S ( 12 , 10 ) =S (10,12) 

DO  1=10,12 
DO  J=10 , 12 
S ( J+3 , 1+3 ) =S ( J , I ) 

S (J+6, 1+6) =S (J, I) 

ENDDO 

ENDDO 

ENDIF 

ENDDO 

DO  I=1,NMAX*6 
DO  K= 1,14 
ESI (K) =0 . DO 
DO  J=1 , 14 

ESI (K) =ESI (K) +B1 ( J, I) *S(J,K) 

ENDDO 

ENDDO 


C 

C*** 


DO  L=I,NMAX*6 
STIFF=0 . DO 
DO  J=1 , 14 

STIFF=STIFF+ESI (J) *B1 (J, L) 
ENDDO 

ES(I,L)=ES(I,L)+STIFF*WT 

ENDDO 


ENDDO 

-k-k-k-k-k-k-k'k-k'k-k'k’k-k-k-k-k-k-k-k-k-k-k'k-k-k’k-k’k-k-k-k-k 
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ENDDO 

C 

ENDDO 

C 

DO  1=1 , 6*NMAX 
DO  J=I+1,6*NMAX 
ES (J, I)=ES (I, J) 
ENDDO 
ENDDO 
C 

RETURN 

END 

C 


APPENDIX  C 

USER  MANUAL  FOR  UFPAC 


UFPAC  is  a finite  element  program  (FORTRAN)  for  linear 
elasto  -static  and  -dynamic  analysis  of  composite  structures. 
The  dynamic  analysis  includes  steady  state,  transient  and 
modal  analyses.  The  initial  stresses  in  the  structure  are 
considered  as  higher  order  effects,  and  the  governing 
eguations  are  linearized  about  the  initial  stress  state.  The 
structural  elements  that  can  be  analyzed  using  the  program  are 
(i)  laminated  beams;  (ii)  laminated  plates;  (iii)  2-D  plane 
solids;  and  (iv)  3-D  solids.  A spring  element  is  also  included 
for  modeling  convenience.  The  laminated  beam  and  plate 
elements  allow  location  of  the  reference  plane  at  an  arbitrary 
distance  from  the  laminate  midplane.  Such  elements  are  called 
offset  elements  and  are  described  in  Chapter  4.  They  provide 
greater  flexibility  in  modeling  laminated  and  sandwich 
structures . 

In  the  static  analysis  option,  displacements,  stresses 
for  continuum  elements  and  stress  resultants  for  structural 
elements,  and  force  resultants  are  calculated  as  specified. 
In  the  modal  analysis  option  modal  loss  factors,  resonant 
frequencies  and  mode  shapes  are  calculated.  In  the  steady 
state  analysis  option,  the  structural  response  for  a forcing 
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function  of  the  type  Felut  is  computed.  The  transient  analysis 
consists  of  finding  the  response  to  a given  set  of  time 
varying  forces.  The  initial  displacements  and  velocities  are 
assumed  to  be  zero.  In  each  case  the  structure  may  be 
prestressed  as  described  in  Chapter  3 . 

The  input  commands  are  read  from  a file  "UFPAC.INP".  The 
laminate  stiffness  matrix,  the  so  called  [ABD]  matrix,  of  the 
beam  and  plate  elements  are  read  from  several  files, 
"ABD1.DAT",  "ABD2.DAT"  etc.  (Appendix  1),  depending  on  the 
problem.  The  data  file  "UFPAC.INP"  expects  commands  in  a 
specific  sequence.  The  output  is  written  onto  a file  called 
"UFPAC* . OUT" . 

The  first  line  of  UFPAC.INP  is  a header,  which  is  printed 
in  the  output  file.  Different  commands  follow  the  header  in  a 
specific  sequence.  Their  names,  definitions  and  data  format 
are  explained  in  the  following  pages.  Commands  are  listed  in 
the  order  they  should  appear  in  the  data  file  UFPAC.INP.  The 
command  FINI  denotes  the  end  of  a set  of  related  group 
commands . 

ANAL ,N1,N2,N3,N4 

ANAL  is  a command  that  specifies  the  type  of  analysis  to  be 
performed.  N1=0:  static  analysis;  Nl=l:  transient  analysis; 
Nl=2:  modal  analysis;  Nl=3:  Direct  Frequency  Response 
Analysis.  N2=0:  No  prestress;  N2=l:  prestressed.  N3=  Number  of 
modes  for  which  frequencies  and  loss  factors  are  to  be 
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calculated  in  the  modal  analysis.  N4=Number  of  master  degrees 
of  freedom  in  the  reduced  modal  analysis.  The  Guyan  Reduction 
Method  (Guyan,  1965)  is  applied  with  automatic  selection  of 
master  degrees  of  freedom.  N4  should  be  set  to  zero  if  no 
reduction  scheme  is  desired.  Maximum  value  for  N4  is  500.  For 
static  and  transient  dynamic  problems  maximum  active  degrees 
of  freedom  allowed  are  3000  with  a half -bandwidth  of  250. 

TYPE ,N1,N2,N3,N4,N5 

TYPE  command  generates  a set  of  elements  of  a specified  type 
with  specified  material  properties.  Nl=Beginning  Element 
Number.  N2=Element  Type.  N3=Material  Number.  N4=Number  of 
elements  to  be  generated.  N5=Element  Number  Increment. 
Example:  TYPE , 49 , 20 , 2 , 9 , 3 

The  above  command  generates  9 elements  49 , 52 , 55 , . . . , 73  . All 
the  nine  elements  are  of  Type  20,  which  is  the  3-D  spring 
element,  and  all  of  them  have  properties  corresponding  to 
Material  2. 

The  following  element  types  are  used  in  the  program. 

N2=10:  Euler-Bernoulli  beam 

N2=20:  Three-Dimensional  uncoupled  springs 

N2=32:  Isoparametric,  3 node,  laminated  beam  element  with 

2 point  integration 

Isoparametric,  3 node,  laminated  beam  element  with 

3 point  integration 


N2=3  3 : 
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N2=62 : 


N2=63 : 


N2=72 : 


N2=73 : 


N2=82 : 


N2=83 : 


N2=92 : 


N2-93 : 


N2=102 : 


N2=103 : 


Isoparametric,  variable  node,  four  sided,  laminated 
plate  element  with  2X2  integration  scheme 
Isoparametric,  variable  node,  four  sided,  laminated 
plate  element  with  3X3  integration  scheme 
Isoparametric,  variable  node,  four  sided, 
orthotropic  plane  element  with  2X2  integration 
scheme 

Isoparametric,  variable  node,  four  sided, 
orthotropic  plane  element  with  3X3  integration 
scheme 

3 node,  7-degree-of-f reedom,  laminated  beam  element 
with  2 point  integration  (quadratic  variation  of  z- 
displacement  and  linear  variation  of  x-displacement 
and  rotation) 

3 node,  7-degree-of-f reedom  laminated  beam  element 
with  3 point  integration  (quadratic  variation  of  z- 
displacement  and  linear  variation  of  x-displacement 
and  rotation) 

Isoparametric,  6 node,  four  sided,  isotropic  plane 
element  with  2X2  integration  scheme 
Isoparametric,  6 node,  four  sided,  isotropic  plane 
element  with  3X3  integration  scheme 
Isoparametric,  variable  node  (from  8 to  20),  six 
faced,  isotropic  solid  element  with  2X2X2 
integration  scheme 

Isoparametric,  variable  node  (from  8 to  20),  six 
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faced,  isotropic  solid  element  with  3X3X3 
integration  scheme 

A FINI  command  is  expected  at  the  end  of  TYPE  commands. 

PROP , N1 , A1 , A2 , A3 , A4 , A5 , A6 

This  command  reads  the  properties  of  material  Nl.  Six 
properties  are  expected.  Some  of  them  can  be  zero.  The 
properties  depend  on  the  type  of  element  the  material  is 
associated  with.  Entries  A1  through  A6  are  expected  in  the 
following  order  depending  upon  the  element  number. 


Element 

type 

10: 

Al-Thickness,  A2-Width,  A3-Young's 

Modulus,  A4-Density,  A5-Loss  Factor 

Element 

type 

20: 

A1-A5  are  the  six  spring  stiffness 
coefficients 

Element 

type 

32 : 

Al-Thickness,  A2-Width,  A3-0ffset, 
Density,  A5-Loss  Factor 

A4- 

Element 

type 

33: 

Same  as  type  32 

Element 

type 

40: 

User  input 

Element 

type 

50: 

User  input 

Element 

type 

62 : 

Al-Density,  A2 -Thickness , A3-Offset, 
Loss  Factor 

A4- 

Element 

type 

63: 

Same  as  type  62 

Element 

type 

72: 

Same  as  type  62 

Element 

type 

73 : 

Same  as  type  62 

Element 

type 

82  : 

Same  as  type  32 

Element 

type 

83  : 

Same  as  type  33 
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Element  type  92: 

Element  type  93: 
Element  type  102: 

Element  type  103: 


Al-Width,  A2-Density , A3-Young's  Modulus 
A4-Poisson  Ratio,  A5-Loss  Factor 
Same  as  type  92 

Al-Young's  Modulus,  A2-Poisson  Ratio  A3- 
Density,  A4-Loss  Factor 
Same  as  Type  102 


A FINI  command  is  expected  at  the  end  of  PROP  commands. 


NODE , N1 , A1 , A2 , A3 , A4 , A5 , A6 , N2 , N3 

The  command  NODE  generates  nodes  on  a straight  line  in  the  x- 
y-z  space.  Node  N1  has  the  coordinates  (A1,A2,A3).  A4,A5  and 
A6  are  increment  in  x,y  and  z coordinates  respectively.  N2  is 
the  total  number  of  nodes  generated  including  the  first  node 
Nl.  N2  has  to  be  at  least  egual  to  1.  When  N2=l  no  node 
generation  is  implied.  N3  is  the  node  number  increment  for 
successive  nodes  generated. 

Example : NODE ,12,1.0,2.0,3.0,4.0,5.0,6.0,3,5 

The  above  example  generates  the  following  nodal  information. 

NODE  x v z_ 

12  1.0  2.0  3.0 

17  5.0  7.0  9.0 

22  9.0  12.0  15.0 


NGEN , Nl , N2 , N3 , A1 , A2 , A3 , N4 , N5 

NGEN  is  used  to  generate  a plane  of  nodes  in  the  x-y-z  space 
from  the  nodal  information  of  a line  of  nodes  already 
generated.  One  can  imagine  that  a line  of  nodes  is  swept  on  a 
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plane  to  generate  new  nodes.  The  parent  line  has  nodes  from  N1 
to  N2  with  a node  number  difference  of  N3.  The  line  is  swept 
N4-1  times.  Each  time  the  node  numbers  of  generated  nodes  are 
incremented  by  N5.  The  x,y,  and  z co-ordinates  of  new  nodes 
are  incremented  by  Al,  A2  and  A3  respectively. 

Example : NGEN ,12,22,5,10.0,20.0,30.0,3,20 

generates  the  following  nodes.  The  co-ordinates  of  nodes  on 
the  parent  line  are  taken  from  the  last  example. 

NODE  x v z_ 


12 

1.0 

2 . 0 

3.0 

17 

5.0 

7.0 

9.0 

22 

9.0 

12 . 0 

15.0 

32 

11.0 

22 . 0 

33.0 

37 

15.0 

27.0 

39.0 

42 

19.0 

32.0 

45.0 

52 

21.0 

42.0 

63.0 

57 

25.0 

47 . 0 

69.0 

62 

29.0 

52 . 0 

75.0 

A FINI  command  is  expected  at  the  end  of  NODE  and  NGEN 
commands . 

ELEM, N1,N2,N3,M1, M2 , M3 , M4 , , MN 

ELEM  command  generates  element  connectivity  information.  The 
starting  element  is  Nl.  Total  number  of  elements  generated  are 
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equal  to  N2  including  the  element  Nl.  Ml ,M2 ,M3 ,M4 , . . .MN  are 
the  global  node  numbers  of  the  first  element  Nl.  The  node 
numbers  of  subsequent  elements  differ  by  N3.  The  element 
(Nl+K)  will  have  the  global  nodes  M1+K*N3,  M2+K*N3, 

M3+K*N3 , . . . ,MN+K*N3 . 

Consider  the  example  of  three-node  beam  elements  given  below. 

Example:  ELEM, 10 , 3 , 5 , 1 , 2 , 3 

will  generate  the  following  three  elements 


Element  No. 

Node-1 

Node-2 

Node-3 

10 

1 

2 

3 

11 

6 

7 

8 

12 

11 

12 

13 

EGEN , Nl , N2 , N3 , Ml , M2 , M3 , M4 , . . . ,MMAX,N4 ,N5 

EGEN  is  similar  to  NGEN.  EGEN  generates  new  sets  of  elements 
by  sweeping  a parent  set  of  elements.  The  parent  set  of 
elements  is  assumed  to  be  generated  already  by  using  EGEN.  The 
parent  set  has  elements  Nl,  N1+N3,  N1+2*N3,  Nl+3 *N3 , . . . , N2 . 
This  set  is  swept  N4-1  number  of  times.  The  corresponding 
element  numbers  are  incremented  by  N5  after  each  sweep.  The 
node  numbers  are  incremented  by  Ml,  M2,  M3,  M4,...,MMAX  after 
each  sweep.  MAX  depends  on  the  element  type.  For  a three-node 
beam  element  MAX=3 , whereas  for  a 20-node  solid  element 
MAX=2 0 . In  the  following  the  beam  elements  generated  in  the 
previous  example  are  considered. 

Example : EGEN ,10,12,1,20,20,20,3,4 
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generates  the  following  elements 


Element  No. 

Node-1 

Node-2 

Node-3 

10 

1 

2 

3 

11 

6 

7 

8 

12 

11 

12 

13 

14 

21 

22 

23 

15 

26 

27 

28 

16 

31 

32 

33 

18 

41 

42 

43 

19 

46 

47 

48 

20 

51 

52 

53 

A FINI  command  is  expected  at  the  end  of  ELEM  and  EGEN 
commands . 

DISP, N1 , N2 , N3 , N4 , N5 , N6 , N7 , N8 , N9 

DISP  command  specifies  the  boundary  condition  codes  (BC-Codes) 
for  nodes  (0=fixed,  Infixed) . Six  BC-Codes  have  to  be 
specified.  The  node  N1  has  the  boundary  condition  codes  N2, 
N3 , N4 , N5 , N6 , and  N7 . The  six  BC-Codes  correspond  to  the 
displacements  in  the  x,  y and  z directions,  and  three 
rotations.  The  sign  conventions  for  rotations  do  not  follow 
the  conventional  vector  notations,  but  follow  the  plate  theory 
conventions,  and  they  are  described  under  element  description. 
The  same  BC-Codes  are  imposed  to  another  N8-1  number  of  nodes. 
The  node  numbers  to  which  the  BC-Codes  are  assigned  are:  Nl, 
N1+N9,  Nl+2  *N9 , Nl+3 *N9 , . . . , N1+ (N8-1) *N9 . 
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Example : DISP, 10,1, 0,1, 0,1, 0,3, 20 

imposes  the  following  boundary  conditions. 

NODE  BC-1  BC-2  BC-3  BC-4  BC-5  BC-6 

10  10  10  10 

30  10  10  10 

50 1 0 1 0 1 0 

A FINI  command  is  expected  after  all  DISP  commands 

FORC ,N1,N2,A1,N3,N4 

FORC  command  prescribes  point  loads  in  the  global  co-ordinates 
at  specified  nodes  in  specified  direction.  N1  is  the  node 
number  at  which  load  is  applied;  N2  is  the  degree  of  freedom 
corresponding  to  the  load;  A1  is  the  magnitude  of  force.  The 
same  force  is  applied  at  N3  number  of  nodes  including  Nl.  The 
node  numbers  at  which  the  load  is  repeated  are  Nl,  N1+N4, 
Nl+2  *N4 , Nl+3  *N4 , . . . ,N1+(N3-1) *N4 . 

Example:  FORC , 12 , 3 , 1000 . 0 , 3 , 15 

will  generate  the  following  loads. 


NODE 

DEGREE  OF  FREEDOM 

MAGNITUDE  OF  FORCE 

12 

3 

1000.0 

27 

3 

1000.0 

42 

3 

1000.0 

OMEG , A1 , A2 , A3 , A4 , A5 , A6 

OMEG  command  is  reguired  only  for  rotating  structures.  Al,  A2 
and  A3  are  respectively  the  x,  y and  z coordinates  of  a point 
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on  the  axis  of  rotation.  A4 , A5  and  A6  are  the  x,  y and  z 
components  of  the  angular  velocity  vector  expressed  in  radians 
per  second. 

A FINI  command  is  expected  after  all  FORC  and  OMEG  commands. 
ELFO , N1 , N2 , N3 

ELFO  provides  the  option  for  computing  the  forces  acting  on  an 
element.  The  forces  are  computed  for  N2  number  of  elements: 
Nl,  N1+N3 , N1+2*N3 , . . . ,N1+(N2-1) *N3 . 

TIME , A1 , A2 , A3 , NSTEP 

TIME  specifies  the  parameters  for  time  marching  algorithm. 
Newmark  method  is  used  for  integration  of  the  equations  of 
motion  in  the  time  domain.  A1  is  the  time  step  size.  A2  and  A3 
are  the  parameters  of  the  Newmark  method.  The  minimum  values 
of  A2  and  A3  are  0.5  and  0.25  respectively.  NSTEP  is  the 
number  of  time  steps  for  which  transient  analysis  is 
performed. 

DEFL, Nl , N2 , N3 

DEFL  command  specifies  the  nodes  at  which  displacement  history 
is  printed.  Nl  is  the  beginning  node.  Total  number  of  nodes 
are  N2 . The  nodes  are:  Nl,  N1+N3 , N1+2*N3,  N1+3*N3,..., 

N1+ (N2-1) *N3 . 

Example:  DEFL, 13, 3, 5 

means  that  the  displacement  history  for  nodes  13,  18  and  23 
will  be  printed. 
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RS LT , N 1 , N2  , N 3 

RSLT  command  specifies  the  elements  for  which  the  stress 
resultants  (or  equivalent  measure  depending  on  the  type  of 
element)  has  to  be  printed.  This  command  is  similar  to  DEFL. 
The  elements  specified  are:  Nl,  N1+N3,  N1+2*N3,  N1+3*N3,..., 
N1+ (N2-1) *N3 . 

PRNT , Nl , N2 , N3 

PRNT  command  specifies  the  time  stations  at  which  results  of 
the  transient  dynamic  analysis  have  to  be  printed.  Nl  is  the 
starting  time  station.  N2  is  the  number  of  stations  at  which 
output  is  provided.  N3  is  the  time  step  increment.  The  time 
stations  at  which  the  results  are  printed  are:  Nl,  N1+N3, 

Nl+2  *N3 , N1+3*N3 , . . . ,N1+(N2-1) *N3 . 

A FINI  is  expected  after  this  command. 

LOAD , Nl , Ml , LI , M2 , L2 , M3 , L3 , . . . 

LOAD  command  specifies  the  nodes  at  which  non-zero  forces  are 
applied,  and  their  directions.  Nl  is  the  number  of  nodes  at 
which  non-zero  forces  are  applied.  Ml,  M2,  M3,...  are  the 

nodes  at  which  forces  are  applied.  LI,  L2 , L3  are  the 
corresponding  degrees  of  freedom.  The  command  expects  to  read 
Nl  number  of  pairs  of  data. 

HIST , Ml , A1 , A2 , A3 , . . . ,AN1 

HIST  command  specifies  the  magnitudes  of  forces  applied  at  Nl 
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number  of  nodes  for  time  station  Ml.  The  number  N1  corresponds 
to  that  used  in  LOAD  command,  and  is  the  number  of  forces 
applied.  Several  lines  are  needed  to  input  the  load  history. 
The  time  stations  Ml  have  to  be  specified  in  ascending  order 
in  the  input  data  file.  The  forces  for  missing  stations  are 
computed  by  linear  interpolation.  The  last  HIST  command  should 
contain  the  forces  for  the  last  time  station  NSTEP  specified 
in  the  TIME  command. 


APPENDIX  D 

EXAMPLE  INPUT  COMMAND  FILES 


Two  input  command  files  prepared  by  using  the  commands 
described  in  Appendix  C are  included  in  this  section.  The  two 
command  files  define  the  problems  in  Sections  8.3  and  8.4. 


COMMAND  FILE  1: 

LOSS  FACTORS  AND  FREQUENCIES  OF  A TAPERED  LAMINATED  BEAM 

ANAL, 2, ,8,500 

TYPE, 1,63, 1,20,1 

TYPE, 21, 63, 3, 20, 1 

TYPE, 41, 103, 2, 20,1 

FINI 

PROP, 1, 1580, 0. 00254 ,-0. 00127 ,0.005 
PROP, 3 ,2800,0.000254,0.000127,0.002 
PROP, 2, 0. 002 1E9, 0.49, 970, 1.0 
FINI 

NODE, 1,, -0.025, ,0.05, 0.00075, ,21, 10 

NODE, 3, ,-0.0125, ,0.05, 0.000375,, 21, 10 

NODE, 5, , , , 0. 05, , , 21, 10 

NODE, 9, ,0.025, ,0.05,-0.00075, ,21,10 

NODE, 7, ,0.0125, ,0.05,-0.000375, ,21,10 

NGEN ,1,209,2,0,0,0.000127,2,1 

FINI 

ELEM, 1,10, 2 0,5, 1,2 1,25, 3, 11, 2 3, 15 
egen, 1,10, 1,4, , , , , , , , ,2,10 
egen ,1,20,1,1,,,,,,,,,2,20 

ELEM, 41, 10, 20, 26, 6, 2, 22, 25, 5, 1, 21, 16,4 , 12,24,15,3,11,23 

EGEN ,41,50,1,4,,,,,,, ,,,,,,,,,,,,,2,10 

FINI 

DISP, 1, 0, 0, 0, 0, 0, 0, 10, 1 
DISP, 13,0,0,0,0,0,0,10,20 
disp, 14, 0,0, 0,0, 0,0, 10, 20 
DISP, 17,0,0,0,0,0,0,10,20 
DISP, 18,0,0,0,0,0,0,10,20 
FINI 
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COMMAND  FILE  2: 


Loss  Factors  Frequencies  and  Mode  Shapes  of  a Rotating  Annulus 

ANAL, 2, ,6,400 

TYPE, 1,63, 1,24,1 

type ,25,103,2,24,1 

type, 49,63,3,24,1 

FINI 

PROP, 1, 158  0, 0. 00254, -0. 00127, 0. 005 
prop ,2,0. 002 ld9 ,0.49,970,1.0 
prop, 3,2800, 0. 000254 ,0.000127,0.002 
FINI 
csys, 2 

node ,1,0.100,-90,0,0,15,0,13,2 
node, 41,0.20,-90,0,0,15,0,13,2 
node ,81,0.30,-90,0,0,15,0,13,2 
node, 121,0.4,-90,0,0,15,0,13,2 
node ,161,0.5,-90,0,0,15,0,13,2 
node, 27 , 0.15,-90,0,0,30,0,7,2 
node ,67,0.25,-90,0,0,30,0,7,2 
node ,107,0.35,-90,0,0,30,0,7,2 
node, 147,0.45,-90,0,0,30,0,7,2 
csys , 0 

ngen, 1,185,2, , ,0.000127,2,1 
FINI 

ELEM, 1,4, 4 0,5, 1,4 1,4 5, 3, 27, 4 3, 29 
egen ,1,4, 1,4, 4, 4, 4, 4, 2, 4, 2, ,6, 4 
egen ,1,24,1,1, , , , , , , , ,2,48 

elem, 25, 4, 4 0,4 6, 6, 2, 42, 45, 5, 1,4 1,3 0,4, 28, 44, 29, 3, 27, 43 

egen ,25,28,1,4,,,,, ,,,2,,2,,2, ,2,,  ,,,,6,4 

FINI 

DISP, 1, 0, 0, 0, 0, 0, 0, 26, 1 
disp, 1,0,1,1,0,1,0,5,40 
disp, 2, 0,1, 1,0, 1,0, 5, 40 
disp, 27,0,1,1,0,1,0,4,40 
disp, 28, 0,1, 1,0, 1,0, 4, 40 
disp, 25, 0,1, 1,0, 1,0, 5, 40 
disp, 26, 0,1, 1,0, 1,0, 5, 40 
disp, 39, 0,1, 1,0, 1,0, 4, 40 
disp, 40, 0,1, 1,0, 1,0, 4, 40 
FINI 
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